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Abstract 

This  paper  surveys  methods  for  simplifying  and  approximating  polygonal  surfaces.  A  polygonal  surface  is  a  piecewise- 
linear  surface  in  3-D  defined  by  a  set  of  polygons;  typically  a  set  of  triangles.  Methods  from  computer  graphics,  com¬ 
puter  vision,  cartography,  computational  geometry,  and  other  fields  are  classified,  summarized,  and  compared  both 
practically  and  theoretically.  The  surface  types  range  from  height  fields  (bivariate  functions),  to  manifolds,  to  non¬ 
manifold  self-intersecting  surfaces.  Piecewise-linear  curve  simplification  is  also  briefly  surveyed. 
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1  Introduction 

The  simplification  of  surfaces  has  become  increasingly 
important  as  it  has  become  possible  in  recent  years  to  cre¬ 
ate  models  of  greater  and  greater  detail.  Detailed  sur¬ 
face  models  are  generated  in  a  number  of  disciplines.  For 
example,  in  computer  vision,  range  data  is  captured  us¬ 
ing  scanners;  in  scientific  visualization,  isosurfaces  are 
extracted  from  volume  data  with  the  “marching  cubes” 
algorithm;  in  remote  sensing,  terrain  data  is  acquired 
from  satellite  photographs;  and  in  computer  graphics  and 
computer-aided  geometric  design,  polygonal  models  are 
generated  by  subdivision  of  curved  parametric  surfaces. 
Each  of  these  techniques  can  easily  generate  surface  mod¬ 
els  consisting  of  millions  of  polygons. 

Simplification  is  useful  in  order  to  make  storage,  trans¬ 
mission,  computation,  and  display  more  efficient.  A  com¬ 
pact  approximation  of  a  shape  can  reduce  disk  and  mem¬ 
ory  requirements  and  can  speed  network  transmission.  It 
can  also  accelerate  a  number  of  computations  involving 
shape  information,  such  as  finite  element  analysis,  colli¬ 
sion  detection,  visibility  testing,  shape  recognition,  and 
display.  Reducing  the  number  of  polygons  in  a  model  can 
make  the  difference  between  slow  display  and  real  time 


display. 

A  variety  of  methods  for  simplifying  curves  and  sur¬ 
faces  have  been  explored  over  the  years.  Work  on  this 
topic  is  spread  among  a  number  of  fields,  making  literature 
search  quite  challenging.  These  fields  include:  cartogra¬ 
phy,  geographic  information  systems  (GIS),  virtual  real¬ 
ity,  computer  vision,  computer  graphics,  scientific  visual¬ 
ization,  computer-aided  geometric  design,  finite  element 
methods,  approximation  theory,  and  computational  geom¬ 
etry. 

Some  prior  surveys  of  related  methods  exist,  notably  a 
bibliography  on  approximation  [45],  a  survey  of  spatial 
data  structures  for  curves  and  surfaces  [106],  and  surveys 
of  triangulation  methods  with  both  theoretical  [6]  and  sci¬ 
entific  visualization  [89]  orientations.  None  of  these  sur¬ 
veys  surface  simplification  in  depth,  however. 

The  present  paper  attempts  to  survey  all  previous  work 
on  surface  simplification  and  place  the  algorithms  in  a  tax¬ 
onomy.  In  this  taxonomy,  we  intermix  algorithms  from 
various  fields,  classifying  algorithms  not  according  to  the 
application  for  which  they  were  designed,  but  according 
to  the  technical  problem  they  solve.  By  doing  so,  we 
find  great  similarities  between  algorithms  from  disparate 
fields.  For  example,  we  find  common  ground  between 
methods  for  representing  terrains  developed  in  cartogra¬ 
phy,  methods  for  approximating  bivariate  functions  devel¬ 
oped  in  computational  geometry  and  approximation  the¬ 
ory,  and  methods  for  approximating  range  data  developed 
in  computer  vision.  This  is  not  too  surprising,  since  these 
are  fundamentally  the  same  technical  problem.  By  calling 
attention  to  these  similarities,  and  to  the  past  duplication 
of  work,  we  hope  to  facilitate  cross-fertilization  between 
disciplines. 

Our  emphasis  is  on  methods  that  take  polygonal  sur¬ 
faces  as  input  and  produce  polygonal  surfaces  as  output, 
although  we  touch  on  curved  parametric  surface  and  vol¬ 
ume  techniques.  Our  polygons  will  typically  be  planar  tri¬ 
angles.  Although  surface  simplification  is  our  primary  in¬ 
terest,  we  also  discuss  curve  simplification,  because  many 
surface  methods  are  simple  generalizations  of  curve  meth¬ 
ods. 
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1.1  Characterizing  Algorithms 

Methods  for  simplifying  curves  and  surfaces  vary  in  their 
generality  and  approach  -  among  surface  methods,  some 
are  limited  to  height  fields,  for  example,  while  others  are 
applicable  to  general  surfaces  in  3-D.  To  systematize  our 
taxonomy,  we  will  classify  methods  according  to  the  prob¬ 
lems  that  they  solve  and  the  algorithms  they  employ.  Be¬ 
low  is  a  list  of  the  primary  characteristics  with  which  we 
will  do  so: 

Problem  Characteristics 

Topology  and  Geometry  of  Input:  For  curves,  the 
input  can  be  a  set  of  points,  a  function  y(x),  a 
planar  curve,  or  a  space  curve.  For  surfaces,  the 
input  can  be  a  set  of  points,  samples  of  a  height 
field  z(x,  y )  in  a  regular  grid  or  at  scattered 
points,  a  manifold 1 ,  a  manifold  with  boundary, 
or  a  set  of  surfaces  with  arbitrary  topology  (e.g. 
a  set  of  intersecting  polygons). 

Other  Attributes  of  Input:  Color,  texture,  and  sur¬ 
face  normals  might  be  provided  in  addition  to 
geometry. 

Domain  of  Output  Vertices:  Vertices  of  the  output 
can  be  restricted  to  be  a  subset  of  the  input 
points,  or  they  can  come  from  the  continuous 
domain. 

Structure  of  Output  Triangulation:  Meshes  can 
be  regular  grids,  they  can  come  from  a  hier¬ 
archical  subdivision  such  as  a  quadtree,  or 
they  can  be  a  general  subdivision  such  as  a 
Delaunay  or  data-dependent  triangulation. 

Approximating  Elements:  The  approximating 
curve  or  surface  elements  can  be  piecewise- 
linear  (polygonal),  quadratic,  cubic,  high 
degree  polynomial,  or  some  other  basis 
function. 

Error  Metric:  The  error  of  the  approximation  is 
typically  measured  and  minimized  with  respect 

1 A  manifold  is  a  surface  for  which  the  infinitesimal  neighborhood  of 
every  point  is  topologically  equivalent  to  a  disk.  In  a  triangulated  mani¬ 
fold,  each  edge  belongs  to  two  triangles.  In  a  triangulated  manifold  with 
boundary,  each  edge  belongs  to  one  or  two  triangles. 


to  Z/>  or  Loo  error2.  Distances  can  be  measured 
in  various  ways,  e.g.,  to  the  closest  point  on  a 
given  polygon,  or  closest  point  on  the  entire  sur¬ 
face. 

Constraints  on  Solution:  One  might  request  the 
most  accurate  approximation  possible  using  a 
given  number  of  elements  (e.g.  line  segments 
or  triangles),  or  one  might  request  the  solu¬ 
tion  using  the  minimum  number  of  elements 
that  satisfies  a  given  error  tolerance.  Some 
algorithms  give  neither  type  of  guarantee,  but 
give  the  user  only  indirect  control  over  the 
speed/quality  tradeoff.  Other  possible  con¬ 
straints  include  limits  on  the  time  or  memory 
available. 

Algorithm  Characteristics 

Speed/Quality  Tradeoff:  Algorithms  that  are  opti¬ 
mal  (minimal  error  and  size)  are  typically  slow, 
while  algorithms  that  generate  lower  quality  or 
less  compact  approximations  can  generally  be 
faster. 

Refinement/Decimation:  Many  algorithms  can 
be  characterized  as  using  either  refinement,  a 
coarse-to-fine  approach  starting  with  a  minimal 
approximation  and  building  up  more  and  more 
accurate  ones,  or  decimation,  a  fine-to-coarse 
approach  starting  with  an  exact  fit,  and  dis¬ 
carding  details  to  create  less  and  less  accurate 
approximations. 

1.2  Background  on  Application  Areas 

The  motivations  for  surface  simplification  differ  from  field 
to  field.  Terminology  differs  as  well. 

2  In  this  paper,  we  use  the  following  error  metrics:  We  define  the  Li 
error  between  two  //-vectors  u  and  v  as  |  |u  —  v|  I2  =  [5Z"=i  (ui  ~  f/)2] 

The  Loo  error,  also  called  the  maximum  error,  is  1 1 u  —  v||oo  = 
max''=1  |  Uj  —  it,- 1 .  We  define  the  squared  error  to  be  the  square  of  the  l,i 
error,  and  the  root  mean  square  or  RMS  error  to  be  the  L 2  error  divided 
by  o/n.  Optimization  with  respect  to  the  Lj  and  L^c  metrics  are  called 
least  squares  and  minimax  optimization,  and  we  call  such  solutions  L,— 
optimal  and  L^— optimal,  respectively. 
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Cartography.  In  cartography,  simplification  is  one 
method  among  many  for  the  “generalization”  of  geo¬ 
graphic  information  [86],  In  that  field,  curve  simplifica¬ 
tion  is  called  “line  generalization”.  It  is  used  to  simplify 
the  representations  of  rivers,  roads,  coastlines,  and  other 
features  when  a  map  with  large  scale  is  produced.  It  is 
needed  for  several  reasons:  to  remove  unnecessary  detail 
for  aesthetic  reasons,  to  save  memory/disk  space,  and  to 
reduce  plotting/display  time.  The  principal  surface  type 
simplified  in  cartography  is,  of  course,  the  terrain.  Map 
production  was  formerly  a  slow,  off-line  activity,  but  it 
is  currently  becoming  more  interactive,  necessitating  the 
development  of  better  simplification  algorithms. 

The  ideal  error  measures  for  cartographic  simplification 
include  considerations  of  geometric  error,  viewer  interest, 
and  data  semantics.  Treatment  of  the  latter  issues  is  be¬ 
yond  the  scope  of  this  study.  The  algorithms  summarized 
here  typically  employ  a  geometric  error  measure  based  on 
Euclidean  distance.  The  problem  is  thus  to  retain  features 
larger  than  some  size  threshold,  typically  determined  by 
the  limits  of  the  viewer’s  perception,  the  resolution  of  the 
display  device,  or  the  available  time  or  memory. 

Computer  Vision.  Range  data  acquired  by  stereo  or 
structured  light  techniques  (e.g.  lasers)  can  easily  produce 
millions  of  data  points.  It  is  desirable  to  simplify  the  sur¬ 
face  models  created  from  this  data  in  order  to  remove  re¬ 
dundancy,  save  space,  and  speed  display  and  recognition 
tasks.  The  acquired  data  is  often  noisy,  so  tolerance  of  and 
smoothing  of  noise  are  important  considerations  here. 

Computer  Graphics.  In  computer  graphics  and  the 
closely  related  fields  of  virtual  reality,  computer-aided  ge¬ 
ometric  design,  and  scientific  visualization,  compact  stor¬ 
age  and  fast  display  of  shape  information  are  vital.  For 
interactive  applications  such  as  military  flight  simulators, 
video  games,  and  computer-aided  design,  real  time  perfor¬ 
mance  is  a  very  important  goal.  For  such  applications,  the 
geometry  can  be  simplified  to  multiple  levels  of  detail,  and 
display  can  switch  or  blend  between  the  appropriate  lev¬ 
els  of  detail  as  a  function  of  the  screen  size  of  each  object 
[13,  52],  This  technique  is  called  multiresolution  model¬ 
ing.  Redisplaying  a  static  scene  from  a  moving  viewpoint 
is  often  called  a  walkthrough.  For  off-line,  more  realistic 


simulations  such  as  special  effects  in  entertainment,  real 
time  is  not  vital,  but  reasonable  speed  and  storage  are  nev¬ 
ertheless  important. 

When  3-D  shape  models  are  transmitted,  compression  is 
very  important.  This  applies  whether  the  channel  has  very 
low  bandwidth  (e.g.  a  modem)  or  higher  bandwidth  (e.g. 
the  Internet  backbone).  The  rapid  growth  of  the  World 
Wide  Web  is  spurring  some  of  the  current  work  in  surface 
simplification. 


Finite  Element  Analysis.  Engineers  use  the  finite  ele¬ 
ment  method  for  structural  analysis  of  bridges,  to  simulate 
the  air  flow  around  airplanes,  and  to  simulate  electromag¬ 
netic  fields,  among  other  applications.  A  preprocess  to 
simulation  is  a  “mesh  generation”  step.  In  2-D  mesh  gen¬ 
eration,  the  domain,  bounded  by  curves,  is  subdivided  into 
triangles  or  quadrilaterals.  In  3-D  mesh  generation,  the  do¬ 
main  is  given  by  boundary  surfaces.  Surface  meshes  of  tri¬ 
angles  or  quadrilaterals  are  first  constructed,  and  then  the 
volume  is  subdivided  into  tetrahedra  or  hexahedra.  The 
criteria  for  a  good  mesh  include  both  geometric  fidelity 
and  considerations  of  the  physical  phenomena  being  simu¬ 
lated  (stress,  flow,  etc).  To  speed  simulation,  it  is  desirable 
to  make  the  mesh  as  coarse  as  possible  while  still  resolving 
the  physical  features  of  interest.  In  3-D  simulations,  sur¬ 
face  details  such  as  bolt  heads  might  be  eliminated,  for  ex¬ 
ample,  before  meshing  the  volume.  This  community  calls 
simplification  “mesh  coarsening”. 


Approximation  Theory  and  Computational  Geometry. 

What  is  called  a  terrain  in  cartography  or  a  height  field  in 
computer  graphics  is  called  a  bivariate  function  or  a  func¬ 
tion  of  two  variables  in  more  theoretical  fields.  The  goal  in 
approximation  theory  is  often  to  characterize  the  error  in 
the  limit  as  the  mesh  becomes  infinitely  fine.  In  compu¬ 
tational  geometry  the  goal  is  typically  to  find  algorithms 
to  generate  approximations  with  optimal  or  near-optimal 
compactness,  error,  or  speed  or  to  prove  bounds  on  these. 
Implementation  of  algorithms  and  low  level  practical  op¬ 
timizations  receive  less  attention. 
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2  Curve  Simplification 

Curve  simplification  has  been  used  in  cartography,  com¬ 
puter  vision,  computer  graphics,  and  a  number  of  other 
fields. 

A  basic  curve  simplification  problem  is  to  take  a  poly- 
gonized  curve  with  n  vertices  (a  chain  of  line  segments  or 
“polyline”)  as  input  and  produce  an  approximating  poly- 
gonized  curve  with  m  vertices  as  output.  A  closely  related 
problem  is  to  take  a  curve  with  n  vertices  and  approximate 
it  within  a  specified  error  tolerance. 

Douglas-Peucker  Algorithm.  The  most  widely  used 
high-quality  curve  simplification  algorithm  is  probably  the 
heuristic  method  commonly  called  the  Douglas-Peucker3 
algorithm.  It  was  independently  invented  by  many  peo¬ 
ple  [99],  [31],  [30,  p.  338],  [5,  p.  92],  [125],  [91,  p. 
176],  [3],  At  each  step,  the  Douglas-Peucker  algorithm  at¬ 
tempts  to  approximate  a  sequence  of  points  by  a  line  seg¬ 
ment  from  the  first  point  to  the  last  point.  The  point  far¬ 
thest  from  this  line  segment  is  found,  and  if  the  distance 
is  below  threshold,  the  approximation  is  accepted,  other¬ 
wise  the  algorithm  is  recursively  applied  to  the  two  sub¬ 
sequences  before  and  after  the  chosen  point.  This  algo¬ 
rithm,  though  not  optimal,  has  generally  been  found  to 
produce  the  highest  subjective-  and  objective-quality  ap¬ 
proximations  when  compared  with  many  other  heuristic 
algorithms  [85,  130],  Its  best  case  time  cost4  is  Q(n),  its 
worst  case  cost  is  0(mn),  and  its  expected  time  cost  is 
about  ©(nlogm).  The  worst  case  behavior  can  be  im¬ 
proved,  with  some  sacrifice  in  the  best  case  behavior,  us¬ 
ing  a  ©(«  log«)  algorithm  employing  convex  hulls  [54], 

A  variant  of  the  Douglas-Peucker  algorithm  described 
by  Ballard  and  Brown  [4,  p.  234]  on  each  iteration  splits 
at  the  point  of  highest  error  along  the  whole  curve,  instead 
of  splitting  recursively.  This  yields  higher  quality  approx¬ 
imations  for  slightly  more  time.  If  this  subdivision  tree  is 

^Pronounced,  and  later  spelled,  due  to  name  change,  “Poiker”. 

4  A  function  is  in  0(f(n))  if  it  is  less  than  or  equal  to  cf(n)  as  n  — y  oo, 
for  some  positive  constant  c.  “O”  is  used  for  upper  bounds. 

A  function  is  in  ©(/(;;))  if  it  is  between  c j  /(;;)  and  C2/©)  as ;;  — *  00, 
for  some  positive  constants  ci,  C2. 

A  function  is  in  £2(/(n))  if  it  is  greater  than  or  equal  to  cf(n)  as  n—>  00, 
for  some  positive  constant  c.  “Q”  is  used  for  lower  bounds. 


saved,  it  is  possible  to  dynamically  build  an  approximation 
for  any  larger  error  tolerance  very  quickly  [18]. 

A  potential  problem  is  that  simplification  can  cause  a 
simple  polygon  to  become  self-intersecting.  This  could  be 
a  problem  in  cartographic  applications. 

Faster  or  Higher  Quality  Algorithms.  There  are  faster 
algorithms  than  Douglas-Peucker,  but  all  of  these  are 
generally  believed  to  have  inferior  quality  [84],  One 
such  algorithm  is  the  trivial  method  of  regular  subsam¬ 
pling  ( known  as  the  “nth-point  algorithm”  in  cartography), 
which  simply  keeps  every  kth  point  of  the  input,  for  some 
k,  discarding  the  rest.  This  algorithm  is  very  fast,  but  will 
sometimes  yield  very  poor  quality  approximations. 

Least  squares  techniques  are  commonly  used  for  curve 
fitting  in  pattern  recognition  and  computer  vision,  but  they 
do  not  appear  to  be  widely  used  for  that  purpose  in  cartog¬ 
raphy. 

Polygonal  Boundary  Reduction.  While  the  Douglas- 
Peucker  algorithm  and  its  variants  are  refinement  algo¬ 
rithms,  curves  can  also  be  simplified  using  decimation 
methods.  Boxer  et  al.  [8]  describe  two  such  algorithms 
for  simplifying  2-D  polygons.  The  first,  due  to  Leu  and 
Chen  [75],  is  a  simple  decimation  algorithm.  It  considers 
boundary  arcs  of  2  and  3  edges.  For  each  arc,  it  computes 
the  maximum  distance  between  the  arc  and  the  chord  con¬ 
necting  its  endpoints.  It  then  selects  an  independent  set  of 
arcs  whose  deviation  is  less  than  some  threshold,  and  re¬ 
places  them  by  their  chords.  The  second  algorithm  is  an 
improvement  of  this  basic  algorithm  which  guarantees  that 
the  approximate  curve  is  always  within  some  bounded  dis¬ 
tance  from  the  original.  They  state  that  the  running  time  of 
the  simple  algorithm  is  ©(«),  while  the  bounded-error  al¬ 
gorithm  requires  0(n  +  r2)  time  where  r  is  the  number  of 
vertices  removed. 

Optimal  Approximations.  Algorithms  for  optimal 
curve  simplification  are  much  less  common  than  heuristic 
methods,  probably  because  they  are  slower  and/or  more 
complicated  to  implement.  In  a  common  form  of  optimal 
curve  simplification,  one  searches  for  the  approximation 
of  a  given  size  with  minimum  error,  according  to  some 
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definition  of  “error”.  Typically  the  output  vertices  are 
restricted  to  be  a  subset  of  the  input  vertices.  A  naive, 
exhaustive  algorithm  would  have  exponential  cost,  since 
the  number  of  subsets  is  exponential,  but  using  dynamic 
programming  and/or  geometric  properties,  the  cost  can  be 
reduced  to  polynomial.  The  Lj-optimal  approximation 
to  a  function  y(x)  can  be  found  in  0(mn 2)  time,  worst 
case,  using  dynamic  programming.  Remarkably,  a  slight 
variation  in  the  error  metric  permits  a  much  faster  algo¬ 
rithm:  the  Loo -optimal  approximation  to  a  function  can 
be  found  in  O(n)  time  [63],  using  visibility  techniques 
(see  also  [123,  91]).  When  the  problem  is  generalized 
from  functions  to  planar  curves,  the  complexity  of  the  best 
Loo -optimal  algorithms  we  know  of  jumps  to  0(n2  log  n) 
[63],  These  methods  use  shortest-path  graph  algorithms 
or  convex  hulls.  For  space  curves  (curves  in  3-D),  there 
are  0(n2  log  m)  Loo -optimal  algorithms  [62]. 

Asymptotic  Approximation.  In  related  work,  McClure 
and  de  Boor  analyzed  the  error  when  approximating 
a  highly  continuous  function  y(x)  using  piecewise- 
polynomials  with  variable  knots  [82,  21],  We  discuss 
only  the  special  case  of  piecewise-linear  approximations. 
They  analyzed  the  asymptotic  behavior  of  the  Lp  error  of 
approximation  in  the  limit  as  m,  the  number  of  vertices 
(knots)  of  the  approximation,  goes  to  infinity.  They 
showed  that  the  asymptotic  Lp  error  with  regular  subsam¬ 
pling  is  proportional  to  m~2,  for  any  p.  The  L/;-optimal 
approximation  has  the  same  asymptotic  behavior,  though 
with  a  smaller  constant.  McClure  showed  that  the  spacing 
of  vertices  in  the  optimal  approximation  is  closely  re¬ 
lated  to  the  function’s  second  derivative.  Specifically,  he 
proved  that  as  m  oo,  the  density  of  vertices  at  each  point 
in  the  optimal  Li  approximation  becomes  proportional  to 
|y"(x)  |2/5.  For  optimal  L0 0  approximations,  the  density  is 
proportional  to  |y"(x)|^2.  Also,  as  m  —>■  oo,  all  intervals 
have  equal  error  in  an  L;,-optimal  approximation. 

The  density  property  and  the  balanced  error  property 
described  above  can  be  used  as  the  basis  for  curve  sim¬ 
plification  algorithms  [82],  Although  adherence  to  nei¬ 
ther  property  guarantees  optimality  for  real  simplification 
problems  with  finite  m,  iterative  balanced  error  methods 
have  been  shown  to  generate  good  approximations  in  prac- 
tice[91,p.  181],  Another  caveat  is  that  many  curves  in  na¬ 
ture  do  not  have  continuous  derivatives,  but  instead  have 


some  fractal  characteristics  [80].  Nevertheless,  these  theo¬ 
retical  results  suggest  the  importance  of  the  second  deriva¬ 
tive,  and  hence  curvature,  in  the  simplification  of  curves 
and  surfaces. 

Summary  of  Curve  Simplification.  The  Douglas- 
Peucker  algorithm  is  probably  the  most  commonly 
used  curve  simplification  algorithm.  Most  implemen¬ 
tations  have  0(mn )  cost,  worst  case,  but  typical  cost 
of  ©(ulogw).  An  optimization  with  worst  case  cost  of 
O(nlogn)  is  available,  however.  Optimal  simplification 
typically  has  quadratic  or  cubic  cost,  making  it  impractical 
for  large  inputs. 

3  Surface  Simplification 

Surfaces  are  more  difficult  to  simplify  than  curves.  In  the 
flight  simulator  field,  lower  level  of  detail  models  have 
traditionally  been  prepared  by  hand  [17].  The  results  can 
be  excellent,  but  the  process  can  take  weeks.  Automatic 
methods  are  preferable  for  large  and  dynamic  databases, 
however. 

If  only  a  single  level  of  detail  is  needed,  then  in  some 
cases,  simplification  can  be  obviated  by  simply  avoiding 
generation  of  redundant  data  in  the  first  place.  In  scientific 
visualization,  for  example,  the  marching  cubes  algorithm 
[90]  is  widely  used.  It  generates  many  tiny  triangles  with¬ 
out  testing  for  coplanarity  between  neighbors.  A  more  so¬ 
phisticated  alternative  is  adaptive  polygonization  that  sub¬ 
divides  finely  only  where  the  surface  is  highly  curved  [7]. 
In  computer  aided  geometric  design,  when  polygonizing 
parametric  surfaces,  rather  than  subdivide  and  polygonize 
a  surface  with  a  regular  (u.  v)  grid,  better  results  are  of¬ 
ten  possible  by  subdividing  adaptively  based  on  curvature 
[10]. 

When  simplification  is  needed  however,  one  of  the  al¬ 
gorithms  summarized  below  can  be  used. 

Taxonomy  of  Surface  Simplification  Algorithms.  We 

categorize  algorithms  at  the  highest  level  according  to  the 
class  of  surfaces  on  which  they  operate: 

•  height  fields  and  parametric  surfaces. 
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•  manifold  surfaces,  and 

•  non-manifold  surfaces. 

Within  each  surface  class  we  often  group  algorithms  ac¬ 
cording  to  whether  they  work  by  refinement  or  decima¬ 
tion.  Within  the  subclasses,  methods  are  generally  listed 
chronologically.  We  have  attempted  to  be  fairly  compre¬ 
hensive,  so  consequently  the  good  methods  are  described 
along  with  the  bad.  As  we  summarize  algorithms,  we 
list  their  computational  complexities  and  quote  empirical 
times,  where  known  (of  course,  hardware,  compilers,  lan¬ 
guages,  and  programming  styles  differ  between  individu¬ 
als,  so  we  must  be  careful  when  judging  based  on  this  in¬ 
formation).  Complexities  are  given  in  terms  of  n,  the  num¬ 
ber  of  vertices  in  the  input,  and  m,  the  number  of  vertices 
in  the  output.  Typically,  m  <$C  n. 

3.1  Height  Fields  and  Parametric  Surfaces 

Height  fields  and  parametric  surfaces  are  the  simplest 
class  of  surfaces.  Within  this  class  of  surfaces,  we  divide 
methods  into  the  following  six  sub-classes:  regular  grid 
methods,  hierarchical  subdivision  methods,  feature  meth¬ 
ods,  refinement  methods,  decimation  methods,  and  opti¬ 
mal  methods. 

Regular  grid  methods  are  the  simplest  techniques,  us¬ 
ing  a  grid  of  samples  equally  and  periodically  spaced  in 
x  and  y.  The  hierarchical  subdivision  methods  are  based 
on  quadtree,  k-d  tree,  and  hierarchical  triangulations  us¬ 
ing  a  divide  and  conquer  strategy.  They  recursively  subdi¬ 
vide  the  surface  into  regions,  constructing  a  tree-structured 
hierarchy.  The  next  four  categories  employ  more  general 
subdivision  and  triangulation  methods,  most  commonly 
Delaunay  triangulation.  Feature  methods  select  a  set  of 
important  “feature”  points  in  one  pass  and  use  them  as 
the  vertex  set  for  triangulation.  Refinement  methods  are 
essentially  generalizations  of  the  Douglas-Peucker  algo¬ 
rithm  from  curves  to  surfaces,  where  intervals  are  replaced 
by  triangles  and  splitting  is  replaced  by  re  triangulating. 
They  start  with  a  minimal  approximation  and  use  multiple 
passes  of  point  selection  and  retriangulation  to  build  up  the 
final  triangulation.  Decimation  methods  use  an  approach 
opposite  that  of  refinement  methods:  they  begin  with  a  tri¬ 
angulation  of  all  of  the  input  points  and  iteratively  delete 


vertices  from  the  triangulation,  gradually  simplifying  the 
approximation.  Refinement  methods  thus  work  top-down, 
while  decimation  methods  work  bottom-up.  The  final  cat¬ 
egory,  “optimal  methods”  are  distinguished  more  for  their 
theoretical  focus  than  for  their  method. 

For  many  height  field  simplification  tasks,  the  input  is  a 
height  field  and  the  output  is  a  general  triangulation,  called 
a  triangulated  irregular  network ,  or  TIN,  in  cartography. 
A  TIN  is  a  mesh  of  triangles  where  height  is  a  function  of 
x  and  y :  II (x,  y).  Examples  of  a  height  field  and  general 
triangulation  are  shown  in  Figures  1  and  2. 


3.1.1  Triangulation 

Most  polygonal  surface  simplification  methods  employ 
triangles  as  their  approximating  elements  when  construct¬ 
ing  a  surface.  For  height  fields  and  parametric  surfaces, 
there  is  a  natural  2-D  parameterization  of  the  surface.  Ba¬ 
sic  triangulation  methods  are  described  in  a  2-D  domain, 
or  in  a  3-D  domain  where  height  z  is  a  function  of  x  and  y. 

In  general,  the  topology  of  the  triangulation  can  be  cho¬ 
sen  using  only  the  xy  projections  of  the  input  points,  or  it 
can  be  chosen  using  the  heights  of  the  input  points  as  well. 
The  latter  approach  is  called  data-dependent  triangulation 
[32], 

The  most  popular  triangulation  method  that  does  not  use 
height  values  is  Delaunay  triangulation;  it  is  a  purely  two- 
dimensional  method.  Delaunay  triangulation  finds  the  tri¬ 
angulation  that  maximizes  the  minimum  angle  of  all  tri¬ 
angles,  among  all  triangulations  of  a  given  point  set.  This 
helps  to  minimize  the  occurrence  of  very  thin  sliver  trian¬ 
gles.  Delaunay  triangulations  have  a  number  of  nice  the¬ 
oretical  properties  that  make  them  very  popular  in  com¬ 
putational  geometry.  In  a  Delaunay  triangulation,  the  cir¬ 
cumscribing  circle  (circumcircle)  of  each  triangle  contains 
no  vertices  in  its  interior  [71],  Delaunay  triangulations 
of  m  points  can  either  be  computed  whole,  using  divide- 
and-conqueror  sweepline  algorithms,  or  incrementally,  by 
inserting  vertices  one  at  a  time,  updating  the  triangula¬ 
tion  after  each  insertion  [48].  The  former  approach  has 
cost  O(mlogm),  while  the  latter,  incremental  method  has 
worst  case  cost  of  0(m2).  Typical  costs  for  the  incremen¬ 
tal  approach  are  much  better  than  quadratic,  however. 
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Figure  1 :  Top  view  of  a  regular  grid  triangu¬ 
lation  of  65  x  65  height  field. 


Figure  2:  A  triangulation  using  512  vertices 
approximating  the  height  field. 


Sometimes  equilateral  triangles  are  not  optimal,  and 
maximization  of  the  minimum  angle  is  not  the  appropriate 
goal.  Triangulation  methods  that  attempt  to  optimize  the 
approximation  of  z  or  other  data  associated  with  the  trian¬ 
gulation  are  called  data-dependent triangulation  methods. 
Several  researchers  have  shown  that  slivers  can  be  good 
when  the  surface  being  approximated  is  highly  curved  in 
one  direction,  but  not  the  other  [102,  88,  20,  32],  Such 
slivers  would  not  be  generated  by  Delaunay  triangulation, 
which  minimizes  slivers  by  tending  to  choose  “fat”  trian¬ 
gles. 

3.1.2  Regular  Grid  Methods 

The  simplest  method  for  approximation  of  surface  grids 
is  regular  subsampling,  in  which  the  points  in  every  Ath 
row  and  column  are  kept  and  formed  into  a  grid,  and  all 
other  points  are  discarded.  Regular  grids  are  also  known 
as  uniform  grids,  and  sometimes  the  term  DEM  (digital  el¬ 
evation  model)  is  used  in  the  specific  sense  of  a  regular 
grid  terrain  model.  As  with  curves,  regular  subsampling 
is  simple  and  fast,  but  low  quality,  since  the  points  dis¬ 
carded  might  be  the  most  important  ones.  The  results  are 


improved  if  a  low  pass  filter  [5 1  ]  is  run  across  the  data  be¬ 
fore  subsampling,  but  this  still  does  not  fix  the  basic  prob¬ 
lem  with  this  method,  its  non-adaptive  nature. 

Kumler  94.  An  extensive  comparison  of  regular  grids 
(DEMs)  and  general  triangulations  (TINs),  and  the 
space/error  tradeoffs  between  them,  was  done  by  Kumler 
[70],  Fie  concluded,  surprisingly,  that  for  a  given  amount 
of  storage  space,  regular  grids  approximate  terrains 
better  than  general  triangulations.  His  comparison  seems 
biased  against  general  triangulations,  however,  since 
he  compares  models  of  equal  memory  size,  not  equal 
rendering  time,  and  the  simplification  algorithms  he  uses 
are  not  the  best  known  [39],  Kumler  assumes  that  general 
triangulations  require  three  to  ten  times  the  memory  of 
regular  grids  with  the  same  number  of  vertices. 

Pyramids.  Regular  subsampling  can  be  done  hierarchi¬ 
cally,  forming  a  pyramid  of  samples  [131,  121],  Despite 
the  wealth  of  research  on  hierarchical  triangulations  and 
TINs,  pyramids  are  probably  the  most  widely  used  type 
of  multiresolution  terrain  model  in  the  simulator  commu- 
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nity  [16,  17]  and  in  the  visualization/animation  commu¬ 
nity  [60],  because  of  their  simplicity  and  compactness. 

3.1.3  Hierarchical  Subdivision  Methods 

Hierarchical  subdivision  methods  construct  a  triangula¬ 
tion  by  recursively  subdividing  a  surface.  They  are  the 
adaptive  form  of  pyramids.  The  hierarchical  pattern  of 
subdivision,  even  if  not  stored  explicitly  in  the  data  struc¬ 
ture,  forms  a  tree,  each  node  of  which  has  no  more  than  one 
parent.  With  Delaunay  triangulation  and  other  general  tri¬ 
angulation  algorithms,  the  topology  is  not  hierarchical,  be¬ 
cause  a  triangle  might  have  multiple  parents.  Hierarchical 
subdivision  methods  are  generally  fast,  simple,  and  they 
facilitate  multiresolution  modeling.  In  perspective  scenes 
where  nearby  portions  of  the  terrain  require  more  detail 
than  distant  regions,  the  hierarchy  facilitates  rendering  at 
adaptive  levels  of  detail.  Nearby  portions  are  drawn  at  at  a 
fine  level,  while  distant  regions  are  drawn  at  a  coarse  level. 
The  penalty  for  their  simplicity  and  speed  is  that  hierar¬ 
chical  subdivision  methods  typically  yield  poorer  quality 
approximations  than  more  general  triangulation  methods. 

Quadtrees  and  k-d  Trees.  Terrains  and  parametric  sur¬ 
faces  are  easily  simplified  using  adaptive  quadtree  and  k-d 
tree  subdivision  methods  [106, 105].  DeHaemer  and  Zyda 
used  quadtree  and  k-d  tree  splitting  at  the  point  of  max¬ 
imum  error  within  each  cell  to  approximate  general  3-D 
surfaces  described  by  a  grid  [27].  Taylor  and  Barrett  de¬ 
scribed  a  similar  method  for  terrains  [122].  Von  Herzen 
and  Barr  discuss  a  method  for  crack-free  adaptive  triangu¬ 
lation  of  parametric  surfaces  using  quadtrees  [  1 28] .  Gross, 
Gatti,  and  Staadt  have  used  wavelets  to  construct  quadtree 
approximations  of  height  fields  [44].  With  an  unoptimized 
implementation,  they  were  able  to  simplify  a  256  x  256  ter¬ 
rain  in  about  2  seconds  on  an  SGI  Indy. 

Gomez-Guzman  79.  A  quaternary  triangulation 
method  for  height  field  approximation  was  proposed  by 
Gomez  and  Guzman  [42].  In  their  method,  each  triangle 
is  recursively  subdivided  into  four  subtriangles  until 
a  maximum  error  tolerance  is  met.  To  subdivide  each 
triangle,  a  “significant”  point  near  the  midpoint  of  each 
edge  is  chosen  (in  some  unspecified  way),  and  the  triangle 


is  split  into  four  nearly  congruent  triangles  (Figure  3). 
Since  the  new  vertices  are  not  constrained  to  lie  on  the 
edges,  however,  the  surface  develops  unsightly  cracks, 
rendering  the  method  unsuitable  for  most  purposes. 

De  Floriani-Falcidieno-Nagy-Pienovi  84.  In  1984,  De 
Floriani  el  al.  published  a  hierarchical  ternary  triangula¬ 
tion  method  in  which  points  are  inserted  in  triangle  inte¬ 
riors  and  each  triangle  is  split  into  three  subtriangles  by 
adding  edges  to  its  vertices  [23],  No  edge  swapping  is 
done  (Figure  4).  Consequently,  all  of  the  initial  edges 
remain  in  the  triangulation  forever,  most  notably  the  di¬ 
agonal  across  the  entire  grid  rectangle,  leading  to  spuri¬ 
ous  knife-edge  ridges  and  valleys  through  the  terrain.  The 
flaws  of  this  method  make  it  unacceptable. 

Schmitt  85.  Schmitt  and  Gholizadeh  simplified  a  grid 
with  rectangular  topology  in  3-D  using  a  triangulated  sur¬ 
face  [112].  Their  method  is  similar  to  that  of  Faugeras 
el  al.  [34],  described  later.  Having  the  input  points  in  a 
grid  allows  the  partition  of  points  into  triangles  to  be  done 
in  a  two  dimensional  parametric  space.  The  method  be¬ 
gins  with  a  small  number  of  triangles  and  repeatedly  splits 
those  triangles  whose  associated  input  points  are  above  the 
error  tolerance.  Triangles  are  subdivided  into  2—4  subtri¬ 
angles  by  splitting  one,  two,  or  three  edges  of  the  triangle. 
Triangle  splitting  is  done  in  no  particular  order.  They  re¬ 
port  that  simplifying  a  grid  of  n  —  288  x  360  points  down 
to  about  m  —  3,  500  points  takes  1.5  hours  on  a  DEC  VAX 
780.  The  computational  complexity  of  their  algorithm  is 
O(mn). 

Scarlatos-Pavlidis  92a.  The  hierarchical  triangulation 
algorithm  for  height  fields  developed  by  Scarlatos  and 
Pavlidis  employs  a  recursive  triangulation  approach  [108, 
107].  Their  method  begins  with  a  minimal  triangulation 
(typically  two  triangles)  as  level  of  detail  0.  Error  toler¬ 
ances  for  each  level  of  detail  in  the  tree  are  specified  by  the 
user.  To  create  level  i  from  level  i—  1,  the  point  of  highest 
error  along  each  triangle  edge  and  in  each  triangle  interior 
is  found,  those  points  with  error  above  the  threshold  for 
level  i  are  taken  as  new  vertices,  and  each  triangle  is  retri¬ 
angulated  using  one  of  five  simple  subdivision  templates 
(Figure  5).  Passes  of  vertex  selection  and  retriangulation 


Figure  5:  Subdivision  templates  for  Scarlatos  and  Pavlidis’  hierarchical  triangulation. 


for  level  i  are  repeated  until  no  more  candidates  for  that 
level  are  found.  All  levels  of  the  hierarchy  are  retained  in 
the  data  structure,  facilitating  adaptive  display  at  any  de¬ 
sired  detail  level.  In  their  analysis,  Scarlatos  and  Pavlidis 
suggest  that  the  cost  of  the  algorithm  is  0(n  log/;).  Our 
analysis  of  their  algorithm  is  that  their  expected  case  is 
O(n\ogm),  but  if  the  hierarchy  is  very  unbalanced,  the 
worst  case  cost  is  0(mn). 


De  Floriani-Puppo  92.  A  similar  method  was  devel¬ 
oped  by  De  Floriani  and  Puppo  [26].  The  triangle  subdivi¬ 
sion  is  more  general,  however.  To  subdivide  a  triangle  for 
a  given  level  in  the  hierarchy,  a  curve  approximation  algo¬ 
rithm  [4]  is  used  to  add  new  vertices  along  the  edges,  then 
additional  points  are  inserted  in  the  interior  of  the  trian¬ 
gle  until  the  error  threshold  is  met  throughout  the  triangle, 
and  the  interior  of  the  triangle  is  retriangulated  using  De¬ 
launay  triangulation.  The  method  appears  to  have  nearly 
identical  flexibility  and  speed  compared  to  Scarlatos  and 
Pavlidis’  method  [108],  but  it  will  probably  yield  slightly 
better  simplification  for  a  given  error  threshold. 


3.1.4  Feature  Methods 

A  simple,  intuitive  approach  to  height  field  simplification 
is  to  make  one  pass  over  the  input  points,  ranking  each  of 
them  using  some  “importance”  measure,  to  select  the  most 
important  points  as  the  vertex  set,  and  construct  a  triangu¬ 
lation  of  these  points.  Typically,  Delaunay  triangulation 
is  used.  Feature  methods  are  quite  popular  in  cartogra¬ 
phy.  Overall,  our  conclusion  is  that  their  quality  relative 
to  many  of  the  other  methods  is  inferior,  so  we  only  sur¬ 
vey  them  briefly  here. 

Important  points,  also  known  as  “features”  or  “criti¬ 
cal  points”  and  the  edges  between  them,  often  known  as 
“break  lines”,  include  such  topographic  features  as  peaks, 
pits,  ridges,  and  valleys.  The  philosophy  of  many  of  the 
feature  approaches  is  that  some  knowledge  about  the  na¬ 
ture  of  terrains  is  essential  for  good  simplification  [129, 
108].  In  a  feature  approach,  the  chosen  features  become 
the  vertex  set,  and  the  chosen  break  lines  (if  any)  become 
edges  in  a  constrained  triangulation  [6],  The  most  com¬ 
monly  used  feature  detectors  are  2x2  and  3x3  lin¬ 
ear  or  nonlinear  filters,  sometimes  followed  by  a  weeding 
process  that  discards  features  that  are  too  close  together. 
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such  as  a  sequence  of  points  along  a  ridge  line.  Such  ap¬ 
proaches  were  employed  by  Peucker-Douglas  and  Chen- 
Guevara  [92,  12],  Some  methods  examine  larger  neigh¬ 
borhoods  of  points  in  an  attempt  to  measure  importance 
more  globally. 

Southard  91.  One  of  the  more  interesting  feature  meth¬ 
ods  is  Southard’s  [120],  He  uses  the  Laplacian  as  a  mea¬ 
sure  of  curvature.  The  rank  of  each  point’s  Laplacian  is 
computed  within  a  moving  window,  analogous  to  a  me¬ 
dian  filter  in  image  processing,  and  all  points  whose  rank 
is  below  some  threshold  are  selected.  This  is  an  im¬ 
provement  over  the  selection  criteria  of  Peucker-Douglas 
and  Chen-Guevara  cited  earlier,  because  it  is  less  sus¬ 
ceptible  to  noise  and  high  frequency  variations,  but  un¬ 
fortunately,  Southard’s  ranking  approach  tends  to  dis¬ 
tribute  points  roughly  uniformly  across  the  domain,  wast¬ 
ing  points  and  leading  to  inferior  approximations,  in  many 
cases.  After  computing  the  Delaunay  triangulation  of  the 
selected  points.  Southard  performs  a  data-dependent  re¬ 
triangulation,  swapping  edges  where  that  would  reduce  the 
sum  of  the  absolute  errors  along  the  edges  in  the  triangu¬ 
lation. 

3.1.5  Refinement  Methods 

Refinement  methods  are  multi-pass  algorithms  that  begin 
with  a  minimal  initial  approximation,  on  each  pass  they 
insert  one  or  more  points  as  vertices  in  the  triangulation, 
and  repeat  until  the  desired  error  is  achieved  or  the  desired 
number  of  vertices  is  used.  For  input  data  in  a  rectangular 
grid,  the  minimal  approximation  is  two  triangles;  for  other 
topologies,  the  initial  approximation  might  be  more  com¬ 
plex.  Incremental  methods  are  typically  used  to  maintain 
the  triangulation  as  refinement  proceeds. 

To  choose  points,  importance  measures  much  like  those 
of  the  feature  methods  can  be  used.  Whereas  feature  meth¬ 
ods  typically  use  importance  measures  that  are  indepen¬ 
dent  of  the  approximation,  in  refinement  algorithms,  the 
importance  of  a  given  point  is  usually  a  measure  of  the 
error  between  it  and  the  approximation.  For  a  height 
field,  the  most  common  metric  for  the  error  is  simply 
the  maximum  absolute  value  of  the  vertical  error,  the  L0 0 
norm.  This  is  the  error  measure  most  closely  related  to  the 


Douglas-Peucker  algorithm. 

Greedy  Insertion.  We  call  refinement  algorithms  that 
insert  the  point(s)  of  highest  error  on  each  pass  greedy  in¬ 
sertion  algorithms,  “greedy”  because  they  make  irrevoca¬ 
ble  decisions  as  they  go  [15],  and  “insertion”  because  on 
each  pass  they  insert  one  or  more  vertices  into  the  trian¬ 
gulation.  Methods  that  insert  a  single  point  in  each  pass 
we  call  sequential  greedy  insertion  and  methods  that  in¬ 
sert  multiple  points  in  parallel  on  each  pass  we  call  paral¬ 
lel  greedy  insertion.  The  words  “sequential”  and  “paral¬ 
lel”  here  refer  to  the  selection  and  re-evaluation  process, 
not  to  the  architecture  of  the  machine.  Many  variations  on 
the  greedy  insertion  algorithm  have  been  explored  over  the 
years;  apparently  the  algorithm  has  been  reinvented  many 
times. 

Fowler-Little  79.  In  1979,  Fowler  and  Little  published 
a  hybrid  algorithm  that  uses  an  initial  pass  of  feature 
selection  using  2x2  filters  to  “seed”  the  triangulation, 
followed  by  multiple  passes  of  parallel  greedy  insertion 
[37].  On  each  of  these  latter  passes,  for  each  triangle,  the 
point  with  highest  error,  or  candidate  point,  is  found,  and 
all  candidate  points  whose  error  is  above  the  requested 
threshold  are  inserted  into  the  triangulation.  (When  the 
point  of  highest  error  falls  on  an  edge,  they  expand  their 
search  for  the  candidate  to  a  sector  of  the  triangle’s  circum- 
circle,  a  quirk  unique  to  their  algorithm.) 

Fowler  and  Little  discussed  two  methods  for  finding 
candidates.  In  their  exhaustive  search  method,  the  error  at 
each  input  point  is  computed  and  tested  against  the  high¬ 
est  error  seen  so  far  for  that  triangle.  In  the  initial  passes 
of  a  greedy  insertion  method,  the  triangles  are  big,  neces¬ 
sitating  the  testing  of  many  points,  but  in  later  passes  the 
triangles  shrink  and  less  testing  per  triangle  is  required. 
As  a  way  to  speed  the  selection  of  candidates,  they  pro¬ 
pose  an  alternative  method  using  hill-climbing,  in  which  a 
test  point  is  initialized  to  the  center  of  the  triangle,  and  it 
repeatedly  steps  to  the  neighboring  input  point  of  highest 
error  until  it  reaches  a  local  maximum,  where  it  becomes 
the  candidate.  This  latter  method  can  be  much  faster,  espe¬ 
cially  for  the  initial  passes,  but  it  would  also  yield  poorer 
quality  approximations  in  many  cases,  because  the  hill 
climbing  might  fail  to  find  the  global  maximum  within  the 


10 


triangle.  Unfortunately,  Fowler  and  Little  did  not  show  a 
comparison  of  the  results  of  the  two  methods,  and  did  not 
analyze  the  speed  of  their  algorithm.  An  approach  similar 
to  Fowler  and  Little’s  was  very  briefly  described  by  Lee 
and  Schachter  [72], 

De  Floriani-Falcidieno-Pienovi  83.  In  1983,  De  Flo- 
riani  el  al.  presented  a  sequential  greedy  insertion  algo¬ 
rithm  [24,  25].  Their  method  is  purer  than  Fowler  and 
Little’s:  it  does  not  seed  the  triangulation  using  feature 
points,  and  it  inserts  a  single  point  on  each  pass,  not  multi¬ 
ple  points.  Consequently,  the  quality  of  its  approximations 
can  be  higher  than  Fowler  and  Little’s.  The  point  inserted 
in  each  pass  is  the  point  of  highest  absolute  error  from  the 
input  point  set.  To  find  this  point  they  apparently  visit  all 
input  points  on  each  pass,  computing  errors.  Their  paper 
says  that  their  algorithm  has  worst  case  cost  of  0(n2),  but 
too  few  details  of  the  algorithm  or  its  data  structures  are 
provided  to  verify  this.  We  will  refer  to  this  paper  and  al¬ 
gorithm  as  “DeFloriani83”. 

De  Floriani  89.  In  later  work,  De  Floriani  published 
an  algorithm  to  build  a  “Delaunay  pyramid”  [22],  a  hier¬ 
archy  of  Delaunay  triangulations,  using  a  variant  of  her 
1983  greedy  insertion  algorithm  to  construct  each  level 
of  the  pyramid.  Her  1989  paper  describes  the  greedy  in¬ 
sertion  algorithm  in  greater  detail  than  her  earlier  papers 
([24,  25]). 

Each  triangle  stores  the  set  of  input  points  it  contains 
and  the  error  of  its  candidate  point.  On  each  pass,  the  set  of 
triangles  is  scanned  to  find  the  candidate  of  highest  error, 
this  point  is  inserted  using  incremental  Delaunay  triangu¬ 
lation,  and  the  candidates  of  all  the  triangles  in  the  mod¬ 
ified  region  are  recomputed.  Recomputing  the  candidate 
of  a  triangle  requires  calculating  the  error  at  each  point  in 
the  triangle’s  point  set. 

De  Floriani  states  that  the  worst  case  time  cost  to  create 
a  complete  pyramid  of  all  n  points  is  0(n2).  We  believe 
that  the  expected  time  cost  of  her  algorithm,  to  select  and 
triangulate  m  points,  is  0(n  log  m  +  m2)  (compare  to  Al¬ 
gorithm  III  in  [40]). 

Because  point  set  traversal  is  used,  rather  than  triangle 
scan  conversion  [36],  this  algorithm  is  not  limited  to  input 


points  in  a  regular  grid,  as  are  most  height  field  approxima¬ 
tion  algorithms.  The  price  of  this  generality  is  speed;  the 
inner  loops  of  a  set  traversal  method  cannot  be  optimized 
as  much  as  those  of  a  scan  conversion  approach. 

Fleller  90.  Heller  explored  a  hybrid  technique  that  he 
called  “adaptive  triangular  mesh  filtering”  [53,  p.  168]. 
This  technique  is  much  like  Fowler  and  Little’s.  The  prin¬ 
cipal  difference  is  that  the  features  are  chosen  not  with 
a  fixed-size  local  filter  but  by  checking  a  variable-sized 
neighborhood  to  determine  if  each  point  is  a  local  ex¬ 
tremum  within  some  height  threshold.  This  feature  selec¬ 
tion  method,  while  more  expensive  than  Fowler  and  Lit¬ 
tle’s,  probably  yields  higher  quality  approximations. 

His  insertion  method  is  sequential,  like  that  of  DeFlo- 
riani83.  He  optimizes  the  algorithm  by  storing  the  set  of 
candidates,  one  candidate  from  each  triangle,  in  a  heap5. 
Below  is  an  excerpt  of  Heller’s  brief  explanation  of  his  al¬ 
gorithm  [53,  p.  168]: 

The  [insertion]  of  a  point  requires  a  local  retrian¬ 
gulation  which  consists  of  swapping  all  neces¬ 
sary  triangles,  and  readjusting  the  [importances] 
of  all  affected  points.  It  is  clear  that  the  time  for 
retriangulation  is  proportional  to  the  number  of 
readjusted  points  and  the  logarithm  of  the  num¬ 
ber  of  queued  points.  It  is,  therefore,  advisable 
to  start  the  process  with  as  many  [feature]  points 
as  possible. 

Due  to  his  optimizations,  Heller’s  algorithm  is  probably 
faster  than  most  others  of  comparable  quality,  such  as 
DeFloriani83,  but  unfortunately,  beyond  the  statements 
quoted  above  he  does  not  analyze  the  speed  of  his  algo¬ 
rithm  theoretically  or  empirically.  It  appears  that  the  ex¬ 
pected  complexity  of  the  greedy  insertion  portion  of  his 
algorithm  is  0((m+n)  log m),  like  Algorithm  III  in  [40], 

Schmitt-Chen  91.  In  order  to  segment  computer  vision 
range  data  into  planar  regions,  Schmitt  and  Chen  use  a 
two  stage  process  called  split-and-merge  [110,  91].  The 

5Christoph  Witzgall  has  also  employed  a  heap.  Personal  communi¬ 
cation.  1994. 
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splitting  stage  is  a  form  of  greedy  insertion  with  Delau¬ 
nay  triangulation  similar  to  DeFloriani83.  The  merging 
stage  joins  together  adjacent  regions  with  similar  normals, 
in  the  process  destroying  the  triangulation,  but  yielding  a 
segmentation  of  the  image.  Their  splitting  stage  approxi¬ 
mated  a  height  field  with  n  =  256 2  points  using  about  m  = 
3,  060  vertices  in  67  seconds  on  a  DEC  VAX  8550. 

Scarlatos-Pavlidis  92a  and  De  Floriani-Puppo  92. 

The  hierarchical  triangulation  methods  of  Scarlatos- 
Pavlidis  [108]  and  De  Floriani-Puppo  [26]  discussed 
earlier  are  analogous  to  greedy  insertion  in  many  ways, 
although  their  triangulations  are  quite  different.  Their 
techniques  will  typically  use  more  triangles  to  achieve  a 
given  error  than  sequential  greedy  insertion  with  Delau¬ 
nay  triangulation,  but  on  the  other  hand,  they  have  the 
advantage  of  a  hierarchy. 

Rippa  92.  Rippa  generalized  the  greedy  insertion  algo¬ 
rithm  of  DeFloriani83  to  explore  data-dependent  triangu¬ 
lation  and  least  squares  fitting  [101]. 

In  place  of  incremental  Delaunay  triangulation,  Rippa’s 
algorithm  computes  a  data-dependent  triangulation  using 
a  version  of  Fawson’s  local  optimization  procedure  [71], 
repeatedly  swapping  edges  around  a  new  vertex  until  the 
global  error  reaches  a  local  minimum.  He  tested  two  defi¬ 
nitions  of  global  error.  The  first  is  a  purely  geometric  mea¬ 
sure:  the  sum  of  the  absolute  values  of  the  angles  between 
normals  of  all  pairs  of  adjacent  triangles  in  the  triangula¬ 
tion,  and  the  second  is  a  simple  L2  measure:  the  sum  of 
squares  of  absolute  vertical  errors  over  all  input  points. 

From  experiments  with  Delaunay  and  data-dependent 
triangulation  on  several  smooth,  synthetic  functions, 
Rippa  concluded  that  data-dependent  triangulation  usu¬ 
ally  yields  more  accurate  approximations  using  a  given 
number  of  vertices  than  Delaunay  triangulation.  The 
angle  criterion  performed  well  in  most  cases,  so  he  mildly 
recommended  it  over  both  the  L2  criterion  and  Delaunay 
triangulation.  Rippa  observed  that  the  L2  criterion  oc¬ 
casionally  allowed  long,  extremely  thin  sliver  triangles 
that  did  not  fit  the  surface  well  to  enter  and  remain  in 
the  triangulation.  The  algorithm  failed  to  eliminate  such 
triangles  because  they  were  so  thin  that  they  contained  no 
input  points,  and  hence  they  contributed  zero  error  to  the 


Z/>  measure. 

The  angle  criterion  also  made  poor  choices  in  some 
cases,  so  Rippa  tried  a  hybrid  scheme  that  on  each  pass 
compares  the  errors  resulting  from  Delaunay  triangula¬ 
tion  and  data-dependent  triangulation  with  the  angle  crite¬ 
rion,  and  updates  using  the  one  with  the  smaller  global  er¬ 
ror.  The  hybrid  scheme  generated  high  quality  approxima¬ 
tions  more  consistently  than  the  other  methods  that  Rippa 
tested.  Unfortunately,  the  hybrid  is  less  elegant,  and  it  ap¬ 
pears  slower  than  the  other  methods.  Margaliot  and  Gots- 
man  reported  an  error  measure  yielding  a  better  fit  than  the 
angle  criterion  [81], 

Rippa  also  explored  least  squares  methods  that  approx¬ 
imate  the  input  points  instead  of  interpolating  them.  The 
(x,  y)  coordinates  of  the  vertices  are  frozen,  but  their 
heights  are  allowed  to  vary,  and  the  combination  of  heights 
that  minimizes  the  global  sum  of  squared  errors  is  found. 
This  involves  solving  a  large,  sparse,  m  x  m  system  of  lin¬ 
ear  equations.  He  found  that  high  quality  results  could  be 
achieved  fairly  efficiently,  on  low-noise  data,  if  the  least- 
squares  fitting  was  done  as  a  post-process  to  greedy  in¬ 
sertion.  His  empirical  tests  on  simple  functions  showed 
that  least  squares  fitting  roughly  halved  the  error  of  the 
standard  interpolative  methods.  Overall,  Rippa’s  methods 
appear  expensive  (data-dependent  triangulation,  particu¬ 
larly  so)  but  the  resulting  approximations  are  higher  qual¬ 
ity  than  those  of  simpler  sequential  greedy  insertion  meth¬ 
ods.  The  least  squares  technique  appears  to  be  particularly 
effective  at  improving  the  approximation. 

Rippa  tested  his  algorithm  on  rather  small  height  fields 
and  did  not  discuss  computational  costs  of  data-dependent 
triangulation  much. 

Polis-McKeown  93.  Polis  and  McKeown  explored 
a  somewhat  parallel  variation  of  the  greedy  insertion 
method  [95].  Their  basic  algorithm,  in  each  pass,  com¬ 
putes  the  absolute  error  at  each  input  point.  The  set  of 
points  of  maximal  absolute  error  is  found,  and  these  are 
inserted  into  the  triangulation,  one  at  a  time,  rejecting 
any  that  are  within  a  tolerance  distance  of  vertices  al¬ 
ready  in  the  triangulation  (see  paper  for  details).  This 
method  might  insert  multiple  points  per  triangle,  unlike 
the  greedy  insertion  algorithms  previously  discussed.  It 
would  typically  insert  fewer  points  per  pass  than  Fowler 
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and  Little’s  algorithm,  however. 

Several  practical  issues  in  the  creation  of  large  ter¬ 
rain  models  for  simulators  are  raised  by  Polis  and  McK- 
eown.  To  facilitate  dynamic  loading  of  the  terrain  as  a 
viewer  roams,  many  display  programs  require  that  ter¬ 
rain  databases  be  broken  into  small  square  blocks  or  “load 
modules”.  This  necessitates  extra  care  along  block  bound¬ 
aries  to  avoid  cracks  between  polygons.  Polis  and  McKe- 
own  also  proposed  selective  fidelity ,  in  which  regions  of 
the  terrain  could  be  assigned  error  weights  according  to 
their  visual  importance,  their  likelihood  of  being  seen,  or 
some  other  criterion.  Thus,  for  example,  for  a  tank  simu¬ 
lator,  one  might  weight  navigable  valleys  more  than  inac¬ 
cessible  mountain  slopes. 

Polis  and  McKeown  tried  a  data-dependent  triangu¬ 
lation  method  involving  summing  the  squares  of  er¬ 
rors  along  all  edges  of  the  triangulation  [94],  much  like 
Southard’s  method.  They  found  Delaunay  triangulation 
to  be  preferable  to  data-dependent  triangulation,  however, 
because  the  former  was  much  faster  [95], 

Polis  and  McKeown’s  algorithm  appears  to  have  an  ex¬ 
pected  cost  of  0(mn )  (like  Algorithm  I  in  [40]).  They  re¬ 
ported  a  compute  time  of  18  hours  to  select  m  —  16,  500 
points  total  from  an  n—  1 , 9792  terrain  broken  into  36  tiles 
on  a  DECstation  5000.  Speed  was  not  the  major  issue  for 
them,  however,  since  they  were  creating  their  TINs  off¬ 
line.  They  later  optimized  their  algorithm  to  select  m  = 
50,  000  points  from  a  terrain  of  n  =  8,  966,  001  points  in 
89  minutes  on  a  DEC  Alpha  [93]. 


Franklin  93.  Franklin  has  released  code  for  a  sequen¬ 
tial  greedy  insertion  algorithm  (PL/I  code  from  1973,  C 
code  from  1993)  [38].  His  algorithmis  quite  similar  to  De- 
Floriani83,  but  optimized  in  a  manner  similar  to  De  Flo- 
riani’s  Delaunay  pyramid  method  ([22]).  With  each  tri¬ 
angle,  Franklin  stores  a  candidate  pointer,  and  he  updates 
only  the  candidates  of  new  or  modified  triangles  on  each 
pass.  He  stores  an  array  of  input  points  with  each  triangle, 
as  in  [22],  so  the  algorithm  is  more  general  but  typically 
slower  than  a  comparable  surface  simplification  algorithm 
limited  to  height  fields. 

Between  his  two  implementations,  Franklin  has  exper¬ 
imented  with  several  triangulation  methods:  swapping  an 


edge  if  it  reduces  the  maximum  error  of  the  approximation, 
swapping  an  edge  if  it  has  shorter  length,  and  Delaunay  tri¬ 
angulation. 

Unfortunately,  Franklin  has  not  published  his  results 
and  conclusions.  By  comparison  to  De  Floriani’s  De¬ 
launay  pyramid  algorithm  and  Algorithm  III  of  [40],  we 
conclude  that  the  expected  cost  of  Franklin’s  algorithm  is 
0(n\ogm  +  nr ).  Franklin’s  program  can  select  m  =  100 
points  from  an  n  = 2572  height  field  in  7  seconds  on  an  SGI 
Indy. 

Puppo-Davis-DeMenthon-Teng  94.  Puppo  et  al.  ex¬ 
plored  terrain  approximation  algorithms  for  the  Connec¬ 
tion  Machine  that  are  parallel  both  in  the  computer  archi¬ 
tecture  sense  and  also  in  the  greedy  insertion  sense  [98]. 
Their  algorithm  is  much  like  that  of  DeFloriani83,  except 
they  insert  all  candidate  points  with  error  above  the  re¬ 
quested  threshold  on  each  pass,  like  Fowler  and  Little. 
They  found  that  the  number  of  points  inserted  on  each  pass 
grew  exponentially,  so  the  number  of  passes  required  to 
insert  m  points  would  typically  be  0(log  in).  On  a  Think¬ 
ing  Machines  CM-2  with  16,384  processors,  they  reported 
compute  times  of  8  seconds  to  select  m  =  379  points  from 
an  n— 1282  terrain  [98],  or  86  seconds  to  select  ni  —  2,  933 
points  from  an  n  —  5122  terrain  [97]. 

The  algorithm  was  parallelized  by  assigning  each  in¬ 
put  point  to  a  different  logical  processor.  Most  of  the  par¬ 
allelization  was  straightforward,  but  parallel  incremental 
triangulation  required  the  use  of  special  mutual  exclusion 
techniques  to  handle  simultaneous  topology  changes  in 
neighboring  triangles. 

Puppo  el  al.  implemented  both  sequential  and  parallel 
greedy  insertion  and  concluded,  surprisingly,  that  the  latter 
is  better.  Our  own  experiments  have  indicated  otherwise 
[40], 

Chen-Schmitt  93.  Chen  and  Schmitt  explored  a  hybrid 
feature/refinement  approach  for  triangulation  of  computer 
vision  range  data  [11].  To  best  approximate  the  step  and 
slope  discontinuities  that  are  common  in  range  data,  they 
first  use  edge  detection  to  identify  significant  discontinu¬ 
ity  features.  These  then  become  constraint  curves  during 
greedy  insertion  of  additional  vertices,  using  either  con- 
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strained  Delaunay  or  data-dependent  triangulation.  Chen 
and  Schmitt  found  that  data-dependent  triangulation  sim¬ 
plified  better  on  surfaces  with  a  preferred  direction,  such 
as  cylinders. 

Silva-Mitchell-Kaufman  95.  A  rather  different  ap¬ 
proach  to  height  field  triangulation  was  proposed  by 
Silva  el  al.  [117].  We  classify  it  here  as  a  refinement 
method,  although  it  is  different  in  spirit  from  the  previous 
methods.  Their  method  uses  greedy  cuts ,  triangulating 
the  domain  from  the  perimeter  inward,  on  each  pass 
“biting”  out  of  the  perimeter  the  triangle  of  largest  area 
that  fits  the  input  data  within  a  specified  maximum  error 
tolerance.  The  method  is  thus  a  generalization  of  greedy 
visibility  techniques  for  curve  simplification  [123,  63], 
and  also  a  form  of  data-dependent  triangulation.  In  a 
comparison  with  Franklin’s  greedy  insertion  algorithm, 
their  unoptimized  program  was  about  two  to  four  times 
slower,  but  produced  triangulations  of  a  given  quality 
using  fewer  vertices.  They  reported  running  times  of 
about  8  seconds  to  select  m=  1, 641  points  from  grids  of 
n  =  1202  points  on  a  one -processor  SGI  Onyx. 

Garland-Heckbert  95.  Our  own  work  in  height  field 
simplification  has  explored  fast  and  accurate  variations  of 
the  greedy  insertion  algorithm  [40,  39], 

We  explored  two  optimizations  of  the  most  basic  greedy 
insertion  algorithm  (as  in  DeFloriani83).  First,  we  ex¬ 
ploited  the  locality  of  mesh  changes,  and  only  recalcu¬ 
lated  the  errors  at  input  points  for  which  the  approxima¬ 
tion  changed,  and  second,  we  used  a  heap  to  permit  the 
point  of  highest  error  to  be  found  more  quickly.  When  ap¬ 
proximating  an  n  point  grid  using  an  m  vertex  triangulated 
mesh,  these  optimizations  sped  up  the  algorithm  from  an 
expected  time  cost  of  0(mn )  to  0((m  +  n)  logm).  We 
were  able  to  approximate  an  n  =  10242  grid  to  high  quality 
using  1%  of  its  points  in  about  21  seconds  on  a  150  MHz 
SGI  Indigo2. 

We  also  explored  a  data-dependent  greedy  insertion 
technique  similar  to  Rippa’s  method.  We  found  an  algo¬ 
rithm  that  yielded,  in  a  fairly  representative  test,  a  solu¬ 
tion  with  88%  the  error  of  Delaunay  greedy  insertion  at  a 
cost  of  about  3-4  times  greater.  Source  code  for  these  al¬ 
gorithms  is  available. 


In  that  paper,  we  propose  several  ideas  for  future  work 
that  could  improve  the  performance  of  the  greedy  inser¬ 
tion  algorithm  in  the  presence  of  cliff  discontinuities,  high 
frequencies,  and  noise. 

Arc/Info  Latticetin.  The  geographic  information  sys¬ 
tem  Arc/Info  sold  by  the  Environmental  Systems  Research 
Institute  (ESRI)  can  approximate  terrain  grids.  Its  “Lat¬ 
ticetin”  command  employs  a  hybrid  feature/refinement  ap¬ 
proach  that  starts  with  a  regular  grid  of  equilateral  trian¬ 
gles  and  refines  it  with  parallel  greedy  insertion  [70,  95], 

3.1.6  Decimation  Methods 

In  contrast  to  refinement  methods,  the  decimation  ap¬ 
proach  to  surface  simplification  starts  with  the  entire  input 
model  and  iteratively  simplifies  it,  deleting  vertices,  trian¬ 
gles,  or  other  geometric  features  on  each  pass.  The  deci¬ 
mation  approach  is  not  so  common  for  height  field  simpli¬ 
fication;  we  will  see  far  more  decimation  methods  in  the 
section  on  manifold  simplification. 

Lee  89.  A  “drop  heuristic”  method  for  simplifying  ter¬ 
rains  was  proposed  by  Lee  [73].  We  call  it  a  vertex  deci¬ 
mation  approach  because  on  each  pass  it  deletes  a  vertex. 
The  algorithm  takes  the  height  field  grid  as  input  and  cre¬ 
ates  an  initial  triangulation  in  which  each  2x2  square  be¬ 
tween  neighboring  input  points  is  split  into  two  triangles 
[73],  On  each  pass,  the  error  at  each  vertex  is  computed 
and  the  vertex  with  lowest  error  is  deleted.  The  error  at  a 
vertex  is  found  by  temporarily  deleting  the  vertex  from  the 
triangulation,  doing  a  local  Delaunay  retriangulation,  and 
measuring  the  vertical  distance  from  the  vertex  to  its  con¬ 
taining  triangle.  The  process  continues  until  the  error  ex¬ 
ceeds  the  desired  level,  or  the  desired  number  of  vertices  is 
reached.  Deletion  in  a  Delaunay  triangulation  can  be  done 
incrementally  to  avoid  excessive  cost  [68]. 

The  drop  heuristic  method  yields  high  quality  approx¬ 
imations,  but  its  computational  cost  and  memory  cost  ap¬ 
pear  very  high.  When  Lee  compared  his  algorithm  to  Chen 
and  Guevara’s  method  and  to  De  Lloriani’s  ternary  trian¬ 
gulation  method  [23],  he  found,  not  surprisingly,  that  his 
method  yielded  superior  results  [74],  The  drop  heuristic 
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method  is  expensive  because  of  the  need  to  visit  each  ver¬ 
tex  on  every  pass.  Its  memory  cost  is  high  because  a  tri¬ 
angulation  with  n  vertices  must  be  created6. 

Scarlatos-Pavlidis  92b.  Scarlatos  and  Pavlidis  explored 
a  method  for  adjusting  a  triangulation  in  order  to  equalize 
the  curvature  of  the  input  data  within  each  triangle  [109], 
extending  McClure’s  and  Pavlidis’  earlier  work  [82,  91, 
83].  Their  algorithm  takes  an  initial  triangulation  and  ap¬ 
plies  three  passes:  shrinking  triangles  with  high  curvature, 
merging  adjacent  coplanar  triangles,  and  swapping  edges 
to  improve  triangle  shape  and  fit.  In  tests,  the  method 
achieved  little  improvement  when  applied  to  the  output 
of  their  hierarchical  triangulation  algorithm  [108,  107]:  in 
most  cases,  the  method  reduced  the  number  of  triangles, 
but  it  also  increased  the  maximum  error  unless  explicit  er¬ 
ror  tests  were  added  [109],  Curvature  equalization  was 
more  successful  at  improving  regular  subsampling  meshes 
[107,  p.  89],  No  unshaded  pictures  of  the  resulting  meshes 
were  given,  however,  so  it  is  difficult  to  compare  the  qual¬ 
ity  of  the  results  to  other  methods. 

Scarlatos  93.  In  addition  to  the  recursive  subdivision 
method  described  earlier,  Scarlatos  also  developed  a  ver¬ 
tex  decimation  method  for  constructing  hierarchical  trian¬ 
gulations  [107].  The  method  begins  with  an  initial  trian¬ 
gulation  and,  to  generate  each  level  of  the  hierarchy,  com¬ 
putes  the  “significance”  of  each  vertex  and  deletes  vertices 
in  increasing  order  of  significance  until  no  more  can  be 
deleted.  Significance  is  an  (unspecified )  function  of  the  er¬ 
ror  between  a  vertex  and  a  weighted  average  of  its  neigh¬ 
bors,  and  the  degree  of  a  vertex.  The  method  is  similar  to 
that  of  Schroeder  el  al. ,  discussed  later,  except  that  Scar¬ 
latos’  method  is  limited  to  height  fields,  and  it  takes  more 
precautions  to  minimize  error  accumulation.  Scarlatos  re¬ 
ported  a  running  time  of  7.75  minutes  to  build  a  complete 
hierarchy  for  about  n  —  5,  900  points  on  a  VAX  8530. 

Hughes-Lastra-Saxe  96.  The  simplification  algorithm 
described  by  Hughes,  Lastra,  and  Saxe  [59]  is  targeted 
towards  simplifying  global  illumination  meshes  resulting 

6 We  find  that  storing  a  triangulation  with  n  vertices  uses  5  to  100 
times  the  memory  of  a  height  field  of  n  points  because  of  the  extra  ad¬ 
jacency  information  required. 


from  radiosity  systems.  Consequently,  the  algorithm  must 
simplify  both  the  mesh  geometry  and  the  color  values  as¬ 
sociated  with  each  mesh  vertex.  They  rejected  a  greedy  in¬ 
sertion  algorithm  because  of  its  inability  to  deal  well  with 
sharp  discontinuities  (i.e.,  shadow  borders).  Instead,  they 
chose  a  combination  of  local  vertex  decimation  and  sim¬ 
plification  envelopes  as  in  [126,  14],  Interestingly,  they 
chose  to  select  vertices  for  removal  at  random  rather  than 
in  order  of  increasing  error.  They  claim  that  this  provides 
more  uniform  meshes,  which  they  believe  to  be  advan¬ 
tageous.  Their  method  also  uses  higher-order  elements 
(quadratic,  cubic,  etc.)  for  reconstructing  the  surface,  a 
possibility  which  most  simplification  methods  ignore. 

3.1.7  Optimal  Methods 

The  error  of  an  optimal  piecewise-linear,  triangulated  ap¬ 
proximation  to  a  smooth  function  of  two  variables  has 
been  analyzed  in  the  limit  as  the  number  of  triangles  goes 
to  infinity.  Nadler  showed  that  the  Li- optimal  approxima¬ 
tion  has  Li  error  proportional  to  m_1  [88]. 

Finding  the  optimal  approximation  of  a  grid  or  surface 
using  triangulations  of  a  subset  of  the  input  points  could  be 
done  by  enumerating  all  possible  subsets  and  all  possible 
triangulations,  but  this  would  take  exponential  time,  and 
it  would  clearly  be  impractical.  As  with  curves,  certain 
problems  in  optimal  surface  approximation  are  well  under¬ 
stood,  while  others  are  not.  It  is  known  that  /-^-optimal 
polygonal  approximation  of  convex  surfaces  is  NP-hard 
(requires  exponential  time,  in  practice)  [19,  9].  This  im¬ 
plies,  of  course,  that  Loo -optimal  approximation  of  height 
fields  and  more  general  surfaces  (in  the  space  of  all  tri¬ 
angulations)  is  also  NP-hard,  since  they  are  a  superset  of 
convex  surfaces.  We  do  not  know  if  there  are  polyno¬ 
mial  time  algorithms  for  optimal  surface  simplification  us¬ 
ing  any  other  error  metric  (such  as  L2),  or  within  a  more 
restricted  class  of  triangulations.  Even  if  some  form  of 
this  problem  permits  an  optimal  algorithm  with  polyno¬ 
mial  time,  it  would  be  surprising  if  it  were  as  fast  as  the 
heuristic  methods  we  have  summarized  above. 

Polynomial  time  algorithms  are  known,  however,  for 
sub-optimal  solutions  with  provable  size  and  quality 
bounds.  If  the  optimal  La 0  solution  for  a  given  error 
tolerance  has  m0  vertices,  there  is  an  0(n7)  algorithm 
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to  find  an  approximation  with  the  same  error  using  m  = 
0(ma  log  m„ )  vertices  [87, 1],  but  this  is  far  too  slow  to  be 
practical  for  large  problems. 

3.2  Manifold  Surfaces 

We  now  turn  our  attention  from  height  fields  and  paramet¬ 
ric  surfaces  to  manifolds  and  manifolds  with  boundary.  In 
general,  the  manifold  can  have  arbitrary  genus  and  be  non- 
orientable7  unless  stated  otherwise.  Manifolds  are  more 
difficult  to  simplify  than  height  fields  or  parametric  sur¬ 
faces  because  there  is  no  natural  2-D  parameterization  of 
the  surface.  Delaunay  triangulation  is  thus  less  easily  ap¬ 
plied.  We  group  manifold  simplification  methods  into  two 
classes:  refinement  methods  and  decimation  methods. 


3.2.1  Refinement  Methods 


Faugeras-Hebert-Mussi-Boissonnat  84.  Faugeras 
et  al.  developed  a  technique  somewhat  similar  to  De 
Floriani’s  1984  algorithm,  but  it  does  not  have  persistent 
long  edges,  and  it  is  applicable  to  the  simplification  of  any 
3-D  triangulated  mesh  of  genus  0,  not  just  height  fields 
[34],  The  method  begins  with  a  pancake-like  two-triangle 
approximation  defined  by  three  vertices  of  the  input  mesh. 
Associated  with  each  triangle  of  the  approximation  is  a 
set  of  input  points.  In  successive  passes,  for  each  triangle 
of  the  approximation,  the  input  point  farthest  from  the 
triangle  is  found,  and  if  the  distance  is  above  threshold, 
the  triangle  is  split  into  3-6  subtriangles  by  inserting 
new  vertices  at  the  interior  point  of  highest  error.  Edges 
common  to  two  subdivided  triangles  are  split  at  their 
points  of  highest  error  (Figure  6).  Splitting  in  this  way 
eliminates  the  long  edges  of  ternary  triangulation. 

During  subdivision,  each  triangle’s  point  set  must  be 
partitioned  into  3-6  subsets.  In  methods  that  are  limited  to 
height  fields,  the  partition  of  input  points  to  subtriangles  is 
done  with  simple  projection  and  linear  splitting.  To  parti¬ 
tion  point  sets  on  a  surface  in  3-D,  Faugeras  el  al.  instead 
split  using  the  shortest  path  along  edges  of  the  input  mesh. 

7  A  manifold  is  orientable  if  its  two  sides  can  be  consistently  labeled 
as  “inside”  and  “outside”.  A  Mobius  strip  is  non-orientable. 


The  method  simplified  an  n  —  2,  000  point  model  in  1 
minute  on  a  Perkin  Elmer  computer.  The  approximations 
generated  were  sometimes  poor,  however,  and  the  method 
had  particular  problems  with  concavities  [96],  A  later  sub¬ 
division  data  structure,  the  “prism  tree”,  addressed  these 
problems  by  recursively  subdividing  surface  points  into 
truncated  pyramidal  volumes  [96]. 

Delingette  94.  A  related  method  for  the  simplification 
of  orientable  manifolds  was  developed  by  Delingette  [28]. 
He  fits  surfaces  to  sets  of  3-D  points  by  minimizing  an  en¬ 
ergy  function  which  is  a  sum  of  an  error  term,  an  edge 
length  term,  and  a  curvature  term.  The  algorithm  starts 
with  a  mesh  that  is  the  dual  to  a  subdivided  icosahedron. 
It  then  iteratively  adjusts  the  geometry,  attempting  to  min¬ 
imize  the  global  energy  [29].  After  a  good  initial  fit  is 
achieved  with  this  fixed  topology,  the  mesh  is  refined.  Re¬ 
gions  of  the  mesh  with  high  curvature,  high  local  fit  er¬ 
ror,  or  elongated  faces  are  subdivided  and  vertices  migrate 
to  points  of  high  curvature  [28].  Delingette  reports  that  it 
takes  2  to  7  minutes  to  approximate  a  set  of  n  —  260,  000 
points  with  a  mesh  of  m  —  1 , 700  vertices  on  a  DEC  Al¬ 
pha.  The  method  is  much  faster  than  the  related  method 
of  Hoppe  el  al.  [58],  but  it  does  not  achieve  comparable 
simplification,  and  it  has  a  number  of  parameters  that  ap¬ 
pear  to  require  careful  tuning. 

Lounsbery-Eck-ef  al.  95.  A  two-stage  method  for  mul¬ 
tiresolution  wavelet  modeling  of  arbitrary  triangulated 
polyhedra  was  developed  by  Lounsbery,  Eck,  et  al.  [76, 
33].  The  method  is  not  limited  to  height  fields  or  even  to 
triangulated  meshes  with  spherical  topology;  it  can  be  ap¬ 
plied  to  any  triangulated  manifold  with  boundary.  The  ap¬ 
proach  first  constructs  a  base  mesh  which  is  a  triangulated 
polyhedron  with  the  same  topology  as  the  input  surface. 
Geodesic-like  distance  measures  are  used  in  this  step,  rem¬ 
iniscent  of  the  method  of  Faugeras  et  al..  It  then  uses  re¬ 
peated  quaternary  subdivision  of  the  base  mesh  to  con¬ 
struct  a  new  mesh  that  approximates  the  input  surface  very 
closely.  A  multiresolution  model  of  the  new  mesh  is  then 
built  using  wavelet  techniques,  after  which  an  approxima¬ 
tion  at  any  desired  error  tolerance  can  be  quickly  gener¬ 
ated.  Eck  et  al.  simplified  a  model  with  about  n  =  35.  000 
vertices  to  m— 5,  400  vertices  in  22  minutes  of  resampling 
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Figure  6:  Subdivision  pattern  of  Faugeras  el  al. . 

plus  5  minutes  of  wavelet  analysis/synthesis,  on  an  SGI 
Onyx.  The  intermediate,  approximating  mesh  had  about 
twice  as  many  vertices  as  the  original. 

While  the  approach  is  very  attractive  for  interactive  sur¬ 
face  design  and  surface  optimization,  it  may  not  be  the  best 
method  for  multiresolution  modeling  of  static  surfaces  be¬ 
cause  of  the  cost  of  resampling.  For  the  approximation  of 
height  fields,  resampling  is  not  needed,  and  simpler  ten¬ 
sor  product  wavelet  techniques  could  be  used  instead  [79] . 

Another  disadvantage  is  that  the  method  does  not  resolve 
creases  at  arbitrary  angles  well,  since  the  final  mesh  sub¬ 
divides  the  triangles  of  the  base  mesh  on  a  regular  grid. 


Figure  7 :  Vertex  decimation.  The  target  vertex  and  its  ad¬ 
jacent  triangles  are  removed.  The  resulting  hole  is  then 
retessellated. 


3.2.2  Decimation  Methods 

The  next  class  of  surface  simplification  algorithms  we  will 
consider  is  decimation  methods:  algorithms  that  start  with 
a  polygonization  (typically  a  triangulation)  and  succes¬ 
sively  simplify  it  until  the  desired  level  of  approximation 
is  achieved.  Most  decimation  algorithms  fall  into  one  of 
the  following  categories: 

vertex  decimation  methods  delete  a  vertex  and  retriangu¬ 
late  its  neighborhood, 

edge  decimation  methods  delete  one  edge  and  two  trian¬ 
gles,  and  merge  two  vertices, 

triangle  decimation  methods  delete  one  triangle  and  three 
edges,  merge  three  vertices,  and  retriangulate  the 
neighborhood,  and 

patch  decimation  methods  delete  several  adjacent  trian¬ 
gles  and  retriangulate  their  boundary. 

Several  variants  of  the  decimation  approach  have  been 
used  for  the  problem  of  simplifying  manifolds,  particu¬ 


larly  for  thinning  the  output  of  isosurface  polygonizers. 

Kalvin  91.  Kalvin  el  al.  developed  a  two  phase  method 
to  create  surface  models  from  medical  data  [65].  The  first 
phase  approximates  a  surface  with  tiny  polygons  using  an 
algorithm  similar  to  marching  cubes  [90],  and  the  second 
phase  then  does  patch  decimation  on  the  model  by  merg¬ 
ing  adjacent  coplanar rectangles.  Since  it  only  merges  pre¬ 
cisely  coplanar  faces,  the  method  does  not  allow  control 
over  the  degree  of  simplification,  so  it  is  quite  limited. 

Schroeder-Zarge-Lorensen  92.  Schroeder  et  al.  devel¬ 
oped  a  general  vertex  decimation  algorithm  primarily  for 
use  in  scientific  visualization  [116].  Their  method  takes 
a  triangulated  surface  as  input,  typically  a  manifold  with 
boundary.  The  algorithm  makes  multiple  passes  over  the 
data  until  the  desired  error  is  achieved.  On  each  pass,  all 
vertices  that  are  not  on  a  boundary  or  crease  that  have  er¬ 
ror  below  the  threshold  are  deleted,  and  their  surrounding 
polygons  are  retriangulated  (see  Figure  7).  The  error  at 


17 


a  vertex  is  the  distance  from  the  point  to  the  approximat¬ 
ing  plane  of  the  surrounding  vertices.  Note  that  errors  are 
measured  with  respect  to  the  previous  approximation,  not 
relative  to  the  input  points,  so  errors  can  accumulate  (this 
flaw  was  fixed  in  later  versions  of  the  algorithm).  Their 
paper  demonstrated  simplifications  of  models  containing 
as  many  as  1,700,000  triangles.  The  computation  time  to 
simplify  a  model  of  n  =  400,  000  vertices  to  m  —  40,  000 
vertices  is  about  14  minutes  on  an  R4000  processor  [115]. 
This  method  uses  significant  memory,  like  Lee’s.  To  con¬ 
serve  memory,  compact  data  structures  were  developed 
[115].  Source  code  for  this  algorithm  is  available  [114]. 

Relative  to  Lee’s  method,  the  technique  of  Schroeder  ef 
al.  is  more  general  since  it  is  not  limited  to  height  fields, 
it  uses  a  less  expensive  and  less  accurate  error  measure, 
and  it  deletes  multiple  vertices  per  pass.  Consequently,  it 
is  faster,  but  probably  has  lower  quality. 

Soucy  and  Laurendeau  92.  To  simplify  manifolds  with 
boundary,  Soucy  and  Laurendeau  also  developed  a  vertex 
decimation  algorithm  [118,  119].  Their  application  was 
the  construction  of  surface  models  from  multiple  range 
views.  On  each  pass,  the  vertex  with  least  error  is  deleted, 
and  its  neighborhood  (the  set  of  adjacent  triangles)  is 
re  triangulated.  The  process  stops  when  the  error  rises 
above  a  specified  tolerance  or  the  desired  size  of  model  is 
achieved. 

To  compute  rigorous  error  bounds,  a  set  of  deleted  ver¬ 
tices  is  stored  with  each  triangle.  We  will  call  these  points 
the  ancestors  of  the  triangle.  To  compute  the  error  at  a 
vertex,  a  temporary  vertex  deletion  and  retriangulation  are 
done.  The  error  of  a  vertex  is  a  measure  of  the  error  that 
would  result  from  its  removal.  More  precisely,  it  is  de¬ 
fined  to  be  the  maximum  distance  between  either  an  an¬ 
cestor  from  the  neighborhood  or  the  vertex  itself  to  the  re¬ 
triangulated  surface.  Deletion  of  a  non-boundary  vertex  is 
considered  legal  if  the  neighborhood  triangles  can  be  pro¬ 
jected  to  2-D  without  foldover. 

To  retriangulate,  Soucy  and  Laurendeau  first  compute 
a  constrained  Delaunay  triangulation  in  a  2-D  projection, 
then  this  triangulation  is  improved  using  a  version  of  Law¬ 
son’s  local  optimization  procedure  [71]  adapted  to  sur¬ 
faces  in  3-D.  To  update  the  data  structures  after  retriangu¬ 
lation,  first  the  ancestor  lists  are  redistributed  among  the 


new  triangles,  then  the  error  of  each  formerly  neighboring 
vertex  is  updated. 

We  can  relate  the  method  to  several  of  its  precursors. 
Like  Lee’s  method,  this  algorithm  does  vertex  decimation 
by  “one  move  lookahead”,  but  unlike  his  technique,  it  is 
not  limited  to  height  fields.  Like  Faugeras  el  al.  and  De 
Floriani  et  al.  (1989),  it  stores  a  point  set  with  each  trian¬ 
gle,  but  unlike  those  methods,  it  is  a  decimation  algorithm, 
and  it  is  more  general:  it  can  simplify  any  manifold  with 
boundary. 

Soucy  and  Laurendeau  estimate  the  expected  complex¬ 
ity  of  their  algorithm  to  be  o(ri  log  (n/(n— m))j .  Their 
method  appears  to  yield  higher  quality  results  than  the 
method  of  Schroeder  et  al..  but  it  is  slower  and  it  uses  more 
memory,  since  it  maintains  lists  of  all  deleted  points.  A  re¬ 
vised  version  of  this  algorithm  is  used  in  the  IMCompress 
software  sold  by  InnovMetric  [64] . 


Turk  92.  Another  method  for  simplifying  a  manifold 
with  boundary  is  due  to  Turk  [124],  This  algorithm  is  not 
a  decimation  method  in  the  same  sense  as  the  previous 
methods,  but  we  list  it  here  because  it  also  starts  with  a  full 
triangulation  and  simplifies. 

Turk’s  algorithm  takes  a  triangulated  surface  as  input, 
sprinkles  a  user-specified  number  of  points  on  these  tri¬ 
angles  at  random,  and  uses  an  iterative  repulsion  proce¬ 
dure  to  spread  the  points  out  nearly  uniformly.  The  points 
remain  on  the  surface  as  they  move  about.  After  these 
points  are  inserted  into  the  original  surface  triangulation, 
the  original  vertices  are  deleted  one  by  one,  yielding  a  tri¬ 
angulation  of  the  new  vertices  that  has  the  same  topology 
as  the  original  surface.  Turk  also  demonstrated  an  im¬ 
proved  variant  of  this  technique  that  groups  points  most 
densely  where  the  surface  is  highly  curved. 

Turk’s  method  appears  to  be  best  for  smooth  surfaces, 
since  it  tends  to  blur  sharp  features8.  Overall,  it  appears 
that  Turk’s  algorithm  is  quite  complex  and  that  it  will  yield 
results  inferior  in  quality  to  the  methods  of  Schroeder  et  al. 
or  Soucy-Laurendeau. 


8William  Schroeder,  SIGGRAPH  ’94  tutorial  talk. 
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Hinker-Hansen  93.  Hinker  and  Hansen  developed  a 
patch  decimation  algorithm  for  use  in  scientific  visualiza¬ 
tion  [55].  It  is  a  one  pass  method  that  first  finds  patches  of 
triangles  with  nearly  parallel  normal  vectors,  and  then  re¬ 
triangulates  each  patch.  The  method  has  O(nlogn)  time 
cost  in  practice.  A  model  with  about  n  —  510,  000  ver¬ 
tices  was  simplified  to  m  =  321, 000  vertices  in  9  min¬ 
utes  on  a  CM-5.  The  method  is  “largely  ineffective  when 
faced  with  surfaces  of  high  curvature”,  however  [55],  It 
appears  to  work  best  on  piecewise-ruled  surfaces:  those 
with  zero  curvature  in  at  least  one  direction,  such  as  cylin¬ 
ders,  cones,  and  planes.  Therefore  the  method  is  not  as 
general  as  that  of  Schroeder  el  al.  or  Soucy-Laurendeau. 

Hoppe-DeRose-Duchamp-McDonald-Stuetzle  93. 

Hoppe  el  al.  developed  an  optimization-based  algorithm 
for  general  3-D  surface  simplification  [58].  Their  method 
takes  a  set  of  points  and  an  initial,  fine  triangulated 
surface  approximation  to  those  points  as  input,  and 
outputs  a  coarser  triangulation  of  the  points  with  the 
same  topology  as  the  input  mesh.  The  method  attempts 
to  minimize  a  global  energy  measure  consisting  of  three 
terms:  a  complexity  term  that  penalizes  meshes  with 
many  vertices,  an  error  term  that  penalizes  geometric 
distance  of  the  surface  from  the  input  points,  and  a  spring 
term  that  penalizes  long  edges  in  the  triangulation.  The 
method  proceeds  in  three  nested  loops,  the  outermost 
one  decreasing  the  spring  constant,  the  middle  one  doing 
an  optimization  over  mesh  topologies,  and  the  inner  one 
doing  an  optimization  over  geometries.  The  topological 
optimization  uses  heuristics  and  random  selection  to 
pick  an  edge  and  either  collapse  it,  split  it,  or  swap  it. 
The  geometric  optimization  uses  nonlinear  optimization 
techniques  to  find  the  vertex  positions  that  minimize  the 
global  error  for  a  given  topology.  Topological  changes 
are  kept  if  they  reduce  the  global  error,  otherwise  they  are 
discarded.  In  other  words,  the  method  makes  repeated 
semi-random  changes  to  the  mesh,  keeping  those  that 
allow  better  fit  and/or  a  simpler  mesh. 

Unlike  most  general  surface  simplification  methods,  the 
method  of  Hoppe  el  al.  does  not  constrain  output  vertices 
to  be  a  subset  of  the  input  points.  Their  method  appears 
to  be  less  sensitive  to  noise  in  the  input  points  than  most 
other  methods  because  of  its  freedom  in  choosing  vertices 
and  because  the  geometric  error  measure  uses  an  Li  norm. 


and  not  an  L0 0  norm. 

Their  method  is  slow,  but  it  is  capable  of  very  good  sim¬ 
plifications.  They  simplified  a  mesh  with  m\  =4,  059  ver¬ 
tices  to  m2 =262  vertices  while  fitting  to  n  — 16,  864  points 
in  46  minutes  on  a  1 -processor  DEC  Alpha.  They  have 
released  their  code.  Their  algorithm  yields  higher  quality 
approximations  than  that  of  Eck  et  al .,  but  it  is  slower  [33]. 

Hamann  94.  A  triangle  decimation  method  was  ex¬ 
plored  by  Hamann  [49].  In  this  algorithm,  triangles  are 
deleted  in  increasing  order  of  weight,  where  weight  is 
the  product  of  “equi-angularity”  and  curvature,  roughly 
speaking.  Thus,  slivers  and  low  curvature  triangles  are 
deleted  first.  The  method  appears  rather  complex,  how¬ 
ever,  since  second  degree  surface  fitting  is  used  to  position 
the  new  vertices,  and  a  number  of  geometric  checks  are  re¬ 
quired  to  prevent  topological  changes. 

Kalvin  94.  In  later  work,  Kalvin  and  Taylor  developed  a 
patch  decimation  method  called  “superfaces”  to  simplify 
manifolds  within  a  given  error  tolerance  [66,  67].  The  al¬ 
gorithm  operates  in  a  single  pass.  This  pass  consists  of 
three  phases.  The  first  phase  segments  the  surface  into  ap¬ 
proximately  planar  patches.  Each  patch  is  found  by  pick¬ 
ing  a  face  at  random  and  merging  in  adjacent  faces  until  the 
patch’s  faces  can  no  longer  be  fit  by  a  plane  within  the  error 
tolerance.  Additional  tests  prevent  degenerate  or  highly 
elongated  patches  from  being  created.  The  second  phase 
simplifies  the  curves  common  to  adjacent  patches  using 
the  Douglas-Peucker  algorithm.  The  third  phase  retrian¬ 
gulates  the  patches  by  subdividing  them  into  star  polygons 
and  then  triangulating  each  star  polygon. 

When  a  face  is  merged  into  a  patch,  the  set  of  feasible 
approximating  planes  ax  +  by  +  cz  +  1  =  0  of  the  patch 
must  be  updated.  This  set  could  be  represented  using  lin¬ 
ear  programming,  as  a  convex  polytope  in  the  3-D  (a,  b ,  c) 
parameter  space  of  planes,  but  the  complexity  of  this  data 
structure  could  grow  quite  large.  Instead,  Kalvin  and  Tay¬ 
lor  use  an  ellipsoidal  approximation  that  supports  constant 
time  updates  and  queries. 

At  a  high  level,  this  method  is  quite  similar  to  Hinker- 
Hansen,  in  that  it  employs  a  single  pass  to  find  nearly 
coplanar  sets  and  then  retriangulates  them.  Hinker- 
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Hansen  define  patches  based  on  angles  between  normal 
vectors,  however,  while  Kalvin-Taylor  define  them  based 
on  distance-to-plane.  Distance  to  plane  is  probably  a  bet¬ 
ter  method  for  defining  patches,  since  it  is  less  sensitive 
to  noise.  Gueziec  reports  that  Kalvin  and  Taylor’s  algo¬ 
rithm  can  simplify  a  model  with  about  n  —  90,  000  vertices 
to  m — 5 ,  000  vertices  in  3  to  5  minutes  on  an  IBM  RS6000. 

Varshney  94.  Using  visibility  techniques  from  compu¬ 
tational  geometry,  Varshney  developed  a  patch  decimation 
algorithm  for  simplifying  orientable  triangulated  mani¬ 
folds  with  boundary  [127,  126],  The  method  has  bounded 
error.  Instead  of  simplifying  in  a  fast,  greedy  manner,  as 
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Figure  8:  A  simple  edge  contraction.  The  highlighted 
edge  is  contracted  into  a  single  point.  The  shaded  trian¬ 
gles  become  degenerate  and  are  removed  during  the  con¬ 
traction. 


most  other  decimation  methods  do,  it  is  much  more  brute 
force,  exhaustively  testing  to  find  the  largest  triangle  to  in¬ 
sert  on  each  pass. 

First,  the  input  surface  is  offset  inwards  and  outwards 
by  a  tolerance  distance  e  to  create  two  offset  surfaces.  All 
triangles  defined  by  three  vertices  of  the  input  surface  are 
checked  for  validity  by  testing  that  they  do  not  intersect 
either  offset  surface  and  that  they  do  not  overlap  previ¬ 
ously  inserted  triangles.  On  each  pass  of  the  algorithm,  the 
valid  triangle  that  “covers”  the  greatest  number  of  previ¬ 
ously  uncovered  input  vertices  is  inserted,  the  old  triangu¬ 
lation  of  this  portion  of  the  surface  is  deleted,  and  small 
triangles  are  added  to  fill  the  cracks  between  the  old  and 
the  new.  The  algorithm  generates  good  approximations 
when  it  works,  but  problems  arise  when  the  offset  surfaces 
collide.  So  far,  the  method  has  not  been  demonstrated  for 


Testing  legality  entails  most  of  the  work  required  to  do 
an  edge  collapse.  The  provisional  new  vertex  is  positioned 
to  fit  the  old  faces  well  and  to  preserve  volume.  During 
simplification,  an  error  radius  is  associated  with  each  ver¬ 
tex.  By  interpolating  spheres  with  these  radii  across  the 
surface,  a  error  volume  is  defined.  At  any  step  during  sim¬ 
plification,  the  error  volume  encloses  the  original  surface. 
When  an  edge  collapse  is  being  considered,  the  error  ra¬ 
dius  for  the  provisional  new  vertex  is  set  so  that  the  new 
error  volume  encloses  the  old  error  volume. 

The  collapse  is  considered  legal  if  it  meets  four  condi¬ 
tions:  (1)  the  topology  of  the  surface  is  preserved,  (2)  the 
normals  of  the  modified  faces  change  little,  (3)  the  new  tri¬ 
angles  are  well  shaped  (not  slivers),  and  (4)  the  error  radius 
for  the  new  vertex  is  below  an  error  threshold. 


simplifications  below  30%  of  the  input  size,  and  it  is  very  Use  of  the  error  volume  could  give  the  user  local  control 
slow.  The  time  costs  of  this  algorithm  and  its  variants  0f  error  tolerance  at  each  vertex.  No  examples  of  this  are 
range  from  <9  Or)  to  0(n6)9  shown  in  the  paper,  however. 


Gueziec  95.  Gueziec  developed  a  method  for  simplify¬ 
ing  orientable  manifolds  that  employs  edge  decimation 
[46],  He  defines  the  edge  collapse ,  or  edge  contraction, 
operator  to  delete  an  edge  and  merge  its  two  endpoints 
into  a  single  vertex  (Figure  8).  Gueziec’s  algorithm  or¬ 
ders  edges  by  “importance”  (in  some  unspecified  way), 
and  makes  a  single  pass  through  the  edges  in  increasing 
order  of  importance,  doing  edge  collapses  where  legal. 

9 Personal  communication,  Pankaj  K.  Agarwal  and  Amitabh  Varsh¬ 
ney,  1995. 


Gueziec  reports  a  time  of  10  minutes  to  simplify  a 
model  with  about  n  —  90,  000  vertices  to  m  —  5,  000  ver¬ 
tices  on  an  IBM  RS6000  model  350.  He  says  that  Kalvin 
and  Taylor’s  algorithm  yields  more  compact  approxima¬ 
tions  for  small  error  tolerances,  but  that  his  algorithm  per¬ 
forms  better  for  large  error  tolerances,  and  that  his  trian¬ 
gles  are  better  shaped.  Closely  related  algorithms  are  il¬ 
lustrated  with  better  pictures  in  another  paper  [47].  In  that 
work,  a  model  with  about  n  =  181 , 000  vertices  was  sim¬ 
plified  to  m  —  26,  000  vertices  in  53  minutes,  on  the  same 
type  of  machine.  The  quality  of  the  resulting  meshes  ap¬ 
pears  good. 
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Gourdon  95.  Gourdon  explored  a  method  for  simplify¬ 
ing  orientable  surface  meshes  resulting  from  surface  re¬ 
construction  [43],  His  algorithm  differs  from  almost  all 
other  simplification  algorithms  in  that  it  does  not  assume 
the  surface  mesh  to  be  a  triangulation.  The  algorithm  is  de¬ 
signed  to  preserve  the  Euler  characteristic10  of  the  model; 
this  implies  that  the  topology  is  preserved.  Topological 
preservation  is  important  for  simplifying  medical  data, 
which  is  the  focus  of  this  technique.  The  algorithm  iter¬ 
atively  removes  edges  based  on  an  unspecified  curvature 
criterion.  Because  the  algorithm  supports  non-triangular 
facets,  no  retessellation  is  required  after  removing  edges. 
Gourdon  observes  that  simply  removing  a  sequence  of 
edges  can  lead  to  undesirable,  irregular  meshes.  To  con¬ 
trol  the  regularity  of  the  tessellation,  he  restricts  the  degree 
of  vertices  to  be  at  most  6  and  facet  may  have  at  most  12 
edges.  Following  simplification,  a  “regularization”  step  is 
performed.  Regularization  attempts  to  improve  the  mesh 
by  moving  points  to  minimize  an  energy  function,  in  this 
case  the  sum  of  squared  edge  lengths.  A  simple  regulariza¬ 
tion  step  would  move  a  vertex  towards  the  barycenter  of  its 
neighbors.  However,  this  can  produce  significant  shrink¬ 
age  of  the  surface.  To  avoid  this,  Gourdon  uses  a  regular¬ 
ization  step  that  moves  the  vertex  towards  the  barycenter, 
but  constrains  the  vertex  to  move  parallel  to  the  average 
plane  of  its  neighbors. 


Klein-Liebich-Strasser  96.  The  algorithm  described  by 
Klein,  Liebich,  and  Strasser  [69]  is  very  similar  to  the 
method  of  Soucy  and  Laurendeau  [119].  It  simplifies  an 
oriented  manifold  by  iteratively  removing  a  vertex  a  retri¬ 
angulating  the  resulting  hole  using  a  constrained  Delaunay 
triangulation.  Each  deleted  vertex  is  linked  to  the  closest 
face  in  the  approximation.  These  links  are  used  to  com¬ 
pute  the  distance  between  the  original  and  approximate 
surfaces.  To  select  a  vertex  for  removal,  each  vertex  is 
tentatively  removed  and  the  additional  error  introduced  by 
the  removal  is  computed.  The  vertex  which  introduces  the 
least  error  is  selected  for  removal.  After  the  vertex  is  re¬ 
moved,  the  links  and  projected  additional  errors  within  its 
neighborhood  must  be  recomputed. 


lnThe  Euler  characteristic  of  a  model  is  defined  as  /  =  F  —  E  +  V 
where  F,  E,  and  V  are,  respectively,  the  number  of  faces,  edges,  and 
vertices. 


Algorri-Schmitt  96.  Algorri  and  Schmitt  developed  an 
algorithm  for  simplifying  closed,  dense  triangulations  re¬ 
sulting  from  surface  reconstruction  [2],  Their  algorithm 
begins  with  a  pre-processing  phase  which  smooths  the  ini¬ 
tial  mesh  by  swapping  edges  based  on  a  G1  -continuity  cri¬ 
terion  as  in  [32].  After  this  initial  smoothing,  every  edge 
whose  dihedral  angle  exceeds  some  user-specified  pla¬ 
narity  threshold  is  classified  as  a  feature  edge.  Each  ver¬ 
tex  is  subsequently  labeled  according  to  its  number  of  inci¬ 
dent  feature  edges.  An  independent  set  of  edges  connect¬ 
ing  “0”  vertices  is  collected,  and  all  the  edges  are  collapsed 
simultaneously.  This  simplification  phase  is  followed  by 
a  smoothing  phase  were  all  non-feature  edges  are  consid¬ 
ered  for  swapping  based  on  the  G1 -continuity  criterion.  If 
further  simplification  is  desired,  edges  are  reclassified  and 
the  process  outlined  above  is  repeated.  Since  only  edges  in 
mostly  planar  regions  are  selected  for  decimation,  the  ba¬ 
sic  step  will  not  simplify  “characteristic  curves”  (e.g.,  the 
edges  of  a  cube)  and  there  will  always  be  a  single  vertex 
left  in  the  midst  of  planar  regions.  Algorri  and  Schmitt  de¬ 
scribe  additional  iterative  steps  which  simplify  these  cases 
separately  from  the  basic  step  outlined  above. 

Ronfard-Rossignac  96.  Another  algorithm  based  on 
edge  collapse  was  described  by  Ronfard  and  Rossignac 
[103].  The  fundamental  observation  underlying  their  al¬ 
gorithm  is  that  each  vertex  in  the  original  model  lies  at  the 
intersection  of  a  set  of  planes,  in  particular,  the  planes  of 
the  faces  that  adjoin  the  vertex.  They  associate  a  set  of 
planes  with  each  vertex;  they  call  this  set  the  zone  of  the 
vertex.  A  vertex’s  zone  is  initialized  to  be  the  set  of  planes 
of  the  adjoining  faces.  The  error  at  a  vertex  is  measured  by 
the  maximum  distance  between  the  vertex  and  the  planes 
in  its  zone.  When  contracting  an  edge,  the  zone  of  the  re¬ 
sulting  vertex  is  the  union  of  the  zones  of  the  original  end¬ 
points.  The  error  of  this  resulting  vertex  characterizes  the 
cost  of  contracting  the  edge.  At  each  iteration,  the  edge  of 
lowest  cost  is  selected  and  contracted.  The  complexity  of 
this  algorithm  would  seem  to  be  0(n  log«). 

Hoppe  96.  The  simplification  algorithm  presented  by 
Hoppe  [56]  for  the  construction  of  progressive  meshes  is 
a  simplified  version  of  the  algorithm  of  Hoppe  el  al.  [58]. 
Rather  than  performing  a  more  general  search,  it  simply 
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Figure  10:  Pair  contraction  joining  unconnected  vertices. 
The  dashed  line  indicates  the  two  vertices  being  contracted 
together. 


Figure  9:  Uniform  vertex  clustering.  Note  the  triangle 
which  as  collapsed  to  a  single  point,  and  the  now  dangling 
edge  at  the  bottom.  Also  note  how  separate  components 
have  been  joined  together. 


selects  a  sequence  of  edge  contractions.  The  algorithm 
uses  essentially  the  same  error  formulation  of  the  earlier 
method,  although  it  is  augmented  to  handle  surface  at¬ 
tributes  such  as  colors.  Hoppe  suggests  that  the  resulting 
meshes  are  just  as  good,  and  perhaps  even  better,  than  the 
results  of  the  more  general  mesh  optimization  algorithm. 

3.3  Non-Manifold  Surfaces 

The  most  general  class  of  surfaces  is  the  non-manifold  sur¬ 
face,  which  permits  three  or  more  triangles  to  share  an 
edge,  and  permits  arbitrary  polygon  intersections.  Rel¬ 
atively  few  surface  simplification  algorithms  can  handle 
models  of  this  generality. 

Rossignac-Borrel  93.  A  very  general  technique  for  sim¬ 
plifying  general  3-D  triangulated  models  was  described 
by  Rossignac  and  Borrel  [104].  They  subdivide  the  ob¬ 
ject’s  bounding  volume  into  a  regular  grid  of  boxes  of  user- 
specified  size.  All  vertices  are  graded  (or  weighted)  ac¬ 
cording  to  some  scheme,  and  all  vertices  within  each  box 
are  merged  together  into  a  new  representative  vertex.  A 
simplified  model  is  then  synthesized  from  these  represen¬ 
tative  vertices  by  forming  triangles  according  to  the  origi¬ 
nal  topology  (see  Figure  9).  This  method  is  extremely  gen¬ 
eral,  as  it  can  operate  on  any  set  of  triangles  (not  just  man¬ 


ifolds),  it  can  achieve  arbitrary  simplification  levels,  and 
it  can  even  eliminate  small  objects  or  otherwise  change 
the  topology  of  a  surface.  Unfortunately,  it  does  not  pre¬ 
serve  detail  well  [52].  When  applied  to  height  fields,  it 
is  roughly  equivalent  to  blurring  followed  by  regular  sub¬ 
sampling.  This  software  is  being  sold  as  part  of  IBM’s  “3D 
Interaction  Accelerator”  [61].  This  method  has  been  ex¬ 
tended  using  octrees  instead  of  regular  grids  [78]. 

Low-Tan  97.  Low  and  Tan  [77]  developed  a  clustering 
algorithm  that  is  intended  to  provide  higher  quality  than 
the  uniform  clustering  described  by  Rossignac  and  Borrel 
while  maintaining  its  generality.  Their  first  improvement 
was  to  suggest  a  better  weighting  criterion.  More  impor¬ 
tantly,  they  replaced  the  uniform  grid  with  a  set  of  cluster 
cells.  These  cells  can  be  any  simple  shape,  such  as  cubes 
or  spheres.  Cells  are  centered  around  their  vertex  of  high¬ 
est  weight.  When  a  vertex  falls  within  the  intersection  of 
multiple  cells,  it  is  placed  in  the  cell  whose  center  is  clos¬ 
est.  In  addition  to  these  algorithmic  improvements,  they 
improved  the  appearance  of  simplified  models  by  render¬ 
ing  stray  edges  as  thick  lines  whose  area  approximates  the 
area  of  the  original  model  in  that  region. 

Garland-Heckbert  97.  We  have  developed  an  algo¬ 
rithm  for  simplifying  surfaces  based  on  iterative  vertex- 
pair  contractions  [41],  A  pair  contraction  is  a  natural  gen¬ 
eralization  of  edge  contraction  (Figure  8)  where  the  vertex 
pair  need  not  be  connected  by  an  edge  (see  Figure  10).  A 
4x4  symmetric  matrix  Q,  is  associated  with  each  vertex 
v,.  The  error  at  the  vertex  is  defined  to  be  vTQv,  and  when 
a  pair  is  contracted,  their  matrices  are  added  together  to 
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form  the  matrix  for  the  resulting  vertex.  We  derive  these 
matrices  to  calculate  the  sum  of  squared  distances  of  the 
vertex  to  a  set  of  planes  (this  is  similar  to  the  error  metric 
of  Ronfard  and  Rossignac  [103]). 

Our  technique  for  tracking  vertex  error  is  quite  efficient, 
and  the  algorithm  is  correspondingly  fast.  The  quality 
of  the  approximations  is  similar  to  those  of  Ronfard  and 
Rossignac,  although  the  algorithm  is  more  general  in  that 
it  can  join  model  components. 

3.4  Related  Techniques 

We  have  focused  on  approximation  of  surfaces  by  poly¬ 
gons,  but  there  has  been  related  work  in  fitting  curved  sur¬ 
faces  to  a  set  of  points  on  a  surface,  and  approximation  of 
volumetric  data.  We  include  a  partial  survey. 

Fitting  a  Curved  Surface  Model.  Polygon  models  for 
curved  surfaces  can  be  bulky.  More  compact  representa¬ 
tions  for  surfaces  are  often  possible  using  curved  surface 
primitives  such  as  piecewise-polynomial  surfaces.  The 
next  class  of  models  beyond  piecewise-linear  surfaces  are 
surfaces  with  tangent  continuity.  Schmitt  and  others  have 
developed  adaptive  refinement  methods  for  fitting  rectan¬ 
gular  Bezier  patches  [113]  and  triangular  Gregory  patches 
[1 1 1]  to  a  grid  of  points  in  3-D.  The  latter  method  is  supe¬ 
rior  to  the  former  because  it  is  better  able  to  adapt  to  fea¬ 
tures  at  an  angle  to  the  grid.  Another  curved  surface  prim¬ 
itive,  the  subdivision  surface,  has  been  fit  to  points  in  3- 
D  by  Hoppe  et  al.,  with  very  nice  results  [57].  Piecewise 
quadratic  surfaces  have  been  fit  to  range  data  using  least 
squares  techniques  [35], 

Fitting  to  a  Volume.  A  generalization  of  the  feature  ap¬ 
proach  to  the  approximation  of  volumes  (scalar  functions 
of  three  variables)  was  explored  by  Hamann  and  Chen 
[50],  They  ranked  points  according  to  an  estimate  of  the 
curvature  of  the  function  fix.  y,  z)  at  each  point,  and  in¬ 
crementally  inserted  vertices  into  a  data-dependent  tetra- 
hedrization,  in  decreasing  order  of  curvature,  until  a  given 
error  tolerance  was  met.  The  errors  for  data-dependent 
tetrahedrization  were  measured  using  Li  or  norms  on 
all  points  inside  each  tetrahedron.  The  surface  decima¬ 


tion  approach  has  also  been  generalized  to  tetrahedriza- 
tions  [100]. 


4  Conclusions 

Surface  simplification  is  not  as  well  understood  as 
curve  simplification.  Whereas  there  appears  to  be  fairly 
widespread  agreement  that  one  algorithm,  Douglas- 
Peucker,  does  a  high  quality  job  of  curve  simplification  at 
acceptable  speeds,  there  is  little  agreement  about  the  best 
approach  for  surface  simplification.  No  thorough  empir¬ 
ical  comparison  of  surface  simplification  methods  has 
been  done  analogous  to  the  studies  for  curves  ([85,  130]). 
Furthermore,  surface  simplification  seems  inherently 
much  more  difficult  than  curve  simplification. 

Why  are  surfaces  so  much  harder?  The  biggest  quali¬ 
tative  difference  we  observe  is  that  curves  inherently  lend 
themselves  to  divide  and  conquer  strategies  like  Douglas- 
Peucker,  since  splitting  a  curve  at  the  point  of  highest  er¬ 
ror  yields  two  curves,  breaking  the  task  into  two  smaller 
subtasks  of  the  same  type.  Splitting  a  surface  at  the  point 
of  highest  error  is  an  ambiguous  concept.  Certain  meth¬ 
ods  arbitrarily  choose  some  way  of  splitting  at  a  point,  as 
with  the  hierarchical  subdivision  methods  that  split  a  tri¬ 
angle  into  three  or  more  subtriangles;  and  other  methods 
abandon  the  divide  and  conquer  strategy  and  employ  the 
more  complex  general  triangulations. 

Our  purpose  has  been  primarily  to  survey  the  existing 
methods,  not  to  evaluate  them,  so  we  offer  few  conclu¬ 
sions  here.  Instead  we  hope  that  this  survey,  by  collecting 
references  and  descriptions  to  the  large  body  of  work  on 
this  topic,  will  draw  attention  to  similar  lines  of  research 
in  disparate  fields,  and  facilitate  future  cross-fertilization. 
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